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1.  INTRODUCTION 


1.1  Motivation  and  Contribution 

Oscillators  are  widely  used  in  RF  circuits  to  generate  reference  signals  [1]. 
Phase  noise,  one  of  the  prime  characteristics  in  the  design  of  RF  and  microwave 
communication  systems  can  be  determined  only  if  the  oscillator  circuit  operates  in 
its  periodic  steady  state.  Therefore,  periodic  steady-state  simulations  of  oscillators 
are  extremely  important.  A  periodic  steady-state  analysis  can  be  performed  using 
the  frequency-domain  harmonic  balance  method  or  the  time-domain  finite- 
difference  or  shooting  method  [2-4]. 

The  harmonic  balance  method  is  a  nonlinear  large-signal  frequency-domain 
technique  that  can  be  used  for  fast  and  accurate  steady-state  analysis  [2-4].  It  can 
achieve  a  superior  accuracy  compared  to  time-domain  methods  due  to  its 
exponential  convergence  behavior  [5].  It  is  also  well  suited  for  RF  system 
simulation  where  the  spectral  and  dynamic  range  requirements  necessitate 
accurate  analysis  of  harmonic  content  with  large  frequency  differences  [6].  Noise 
in  an  oscillator  is  frequently  measured  and  described  in  the  frequency  domain. 
Although  the  harmonic  balance  method  is  well  suited  for  these  applications,  it  is 
susceptible  to  convergence  problems  when  applied  to  oscillator  circuits  [2,  3]. 

The  standard  technique  for  frequency-domain  oscillator  simulation  using  the 
harmonic  balance  method  is  to  apply  a  voltage  probe  which  converts  the  oscillator 
circuit  into  a  set  of  closely  related  numerically  more  efficient  forced  circuits  [7]. 
The  voltage  probe  is  usually  implemented  with  a  two-level  Newton  method.  The 
success  of  the  two-level  method,  however,  depends  on  1)  a  close  initial  guess  for 
the  oscillation  frequency  and  the  probe  voltage,  and  2)  the  convergence  of  the 
bottom-level  probe-forced  circuit  equations.  Depending  on  the  type  of  oscillator 
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the  two-level  Newton  method  fails  to  converge  because  a  good  initial  guess  can 
not  be  provided  (high-Q  oscillators)  [8,  9]  or  the  bottom-level  equations  fail  to 
converge  (highly-nonlinear  oscillators,  e.g.,  ring  oscillators).  For  this  reason, 
alternate  approaches  have  to  be  studied.  In  this  work,  new  algorithms  and 
techniques  have  been  developed  for  the  efficient  and  robust  simulation  of  both 
high-Q  LC  oscillators  and  highly  nonlinear  ring  oscillators. 

A  detailed  analysis  in  [8,  9]  shows  that  a  continuation  method  in 
conjunction  with  the  voltage  probe  of  [7]  can  be  successfully  used  to  simulate  a 
high-Q  oscillator.  A  multistage  continuation  approach  for  oscillator  simulation 
using  a  Hopf  bifurcation  is  also  proposed  in  [10].  However,  the  specific  details  of 
the  continuation  methods  are  not  provided  [8-10].  This  thesis  explores  the  use  of 
homotopy  methods  [11,  12]  for  the  steady-state  simulation  of  high-Q  oscillators  in 
the  frequency  domain. 

Homotopy  methods  are  globally  convergent  numerical  techniques  for  solving 
nonlinear  algebraic  equations.  They  have  been  used  extensively  for  finding  the  dc 
operating  points  of  “difficult-to-solve”  nonlinear  circuits  [12-15].  Recently, 
homotopy  methods  have  also  been  applied  to  the  steady-state  analysis  of 
oscillators  in  the  time  domain  [16-18].  However,  the  use  of  homotopy  methods 
with  frequency-domain  methods  has  not  been  previously  studied.  This  work 
examines  artificial  parameter  homotopy  methods  for  oscillator  steady-state 
simulation  with  the  harmonic  balance  method  [19].  Different  embedding 
techniques  are  compared  and  evaluated.  A  robust  and  efficient  embedding  is 
identified  and  demonstrated  to  work  well  for  a  wide  variety  of  LC  oscillator 
circuits. 

Previous  work  using  the  harmonic  balance  method  for  oscillator  simulation 
has  focused  on  LC  oscillators  which  have  near  sinusoidal  waveforms  [7-9].  Some 
of  these  techniques  can  be  applied  to  simulate  ring  oscillators,  but  due  to  the 
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highly  nonlinear  nature  of  the  circuit,  convergence  problems  exist  [20]. 
Commercial  harmonic  balance  tools  use  results  from  a  time-domain  transient 
analysis  for  simulating  ring  oscillators.  The  initial  guesses  thus  achieved  are 
nearly  the  same  as  the  solution  from  the  frequency-domain  harmonic  balance 
method  and  the  convergence  problem  in  the  frequency-domain  simulation  of  ring 
oscillators  is  not  actually  addressed.  More  importantly  this  approach  cannot  be 
extended  to  the  simulation  of  voltage  controlled  ring  oscillators  in  a  RF  system 
where  transient  analysis  is  no  longer  efficient  due  to  the  presence  of  mixer  circuits. 

In  this  work,  convergence  problems  associated  with  the  frequency-domain 
simulation  of  highly  nonlinear  ring  oscillators  have  been  examined.  Two  efficient 
and  robust  methods,  a  single-delay  cell  method  [21,  22]  and  a  multiple-probe 
method,  are  proposed  for  robust  frequency-domain  ring  oscillator  simulation 
using  the  harmonic  balance  method.  Both  approaches  exploit  the  periodic 
structure  of  a  ring  oscillator  and  have  been  tested  on  a  variety  of  single-ended  and 
differential  ring  oscillators.  They  converge  robustly  for  circuits  where  the 
conventional  methods  fail. 

The  application  of  the  harmonic  balance  method  to  relaxation  oscillator 
simulation  is  also  investigated.  It  is  shown  that  due  to  accuracy  as  well  as 
convergence  problems,  the  standard  harmonic  balance  method  is  not  suitable  for 
relaxation  type  oscillators  which  have  sharp  waveform  transitions. 


1.2  Dissertation  Outline 

This  dissertation  is  organized  as  follows.  Chapter  1  provides  the  motivation 
and  contributions  of  this  work.  In  Chapter  2,  the  frequency-domain  harmonic 
balance  method  is  reviewed.  The  equation  formulation,  solution  strategy,  spectral 
accuracy  and  the  error  mechanisms  for  the  harmonic  balance  method  are 
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introduced.  Oscillator  simulation  using  the  harmonic  balance  method  is  presented 
in  Chapter  3.  First,  oscillator  structures  and  analysis  are  introduced.  Then  the 
application  of  the  harmonic  balance  method  to  simulation  of  oscillators  is 
discussed.  The  concept  of  the  probe  is  introduced  and  its  implementation  by  the 
two-level  method  is  presented.  Convergence  problems  and  challenges  are  then 
described. 

Frequency-domain  simulation  of  high-Q  oscillators  is  presented  in  Chapter 
4.  The  challenges  in  simulating  high-Q  oscillators  are  first  discussed  and  the 
theory  of  homotopy  methods  is  reviewed.  This  is  followed  with  an  algorithm  in 
which  the  homotopy  method  is  applied  to  robust  frequency-domain  high-Q 
oscillator  simulation. 

In  Chapter  5  the  frequency-domain  simulation  of  highly  nonlinear  ring 
oscillators  is  discussed.  The  problems  associated  with  ring  oscillator  simulation 
are  first  introduced.  To  overcome  the  convergence  problems  in  ring  oscillator 
simulation  two  new  methods,  a  single-delay  cell  method  and  a  multiple-probe 
method  are  proposed.  In  the  single-delay  cell  method,  an  equivalent  circuit  with 
only  one  delay  cell  is  used  for  ring  oscillator  simulation.  A  modified  voltage 
probe  technique  is  used  for  convergence.  Implementation  details  are  given,  along 
with  the  limitations  of  this  method.  In  the  multiple-probe  method,  multiple 
probes  are  included  into  the  ring  oscillator  circuit.  A  discussion  on  initial  guesses 
and  computational  effort  is  provided. 

The  application  of  the  harmonic  balance  method  to  simulation  of  relaxation 
oscillators  is  briefly  investigated  in  Chapter  6. 

Chapter  7  presents  the  simulation  results  for  the  new  methods  proposed  in 
this  work.  These  include  simulations  of  high-Q  oscillators  with  the  new 
homotopy-based  algorithm  and  simulation  of  ring  oscillators  using  the  single- 
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delay  cell  and  multiple-probe  methods.  Conclusions  and  future  work  are 
summarized  in  Chapter  8. 
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2.  FREQUENCY-DOMAIN  HARMONIC  BALANCE  METHOD 


2.1  Introduction 

With  the  rapid  growth  of  the  radio  frequency  (RF)  integrated  circuit  market, 
there  is  an  increasing  demand  to  accurately  simulate  system  performance  before 
fabrication  [2-4].  Many  system  characteristics,  such  as  distortion,  power, 
frequency,  and  noise,  can  be  properly  determined  only  if  the  circuit  operates  in  its 
periodic  steady-state  domain.  For  this  reason,  periodic  steady-state  simulations  of 
RF  circuits  are  extremely  important. 

A  trivial  method  for  determining  the  steady- state  response  is  to  run  a 
transient  analysis  and  wait  for  all  the  waveforms  to  settle  to  their  periodic  steady 
state.  However,  this  may  be  inefficient  and,  in  some  cases,  not  suitable  for 
simulating  RF  systems  with  high-Q  filters  and  wide  frequency  separations.  In 
these  cases,  a  long  transient  needs  to  be  simulated  to  make  sure  the  effects  of  the 
initial  conditions  die  out  and  the  time-step  used  by  the  numerical  integration 
algorithm  needs  to  be  chosen  small  enough  to  capture  the  highest  frequency  of 
interest.  Therefore,  a  large  number  of  time  points  have  to  be  simulated.  To 
overcome  the  difficulty  with  the  conventional  transient  simulation,  steady-state 
solution  is  computed  directly  in  the  time  domain  by  the  finite-difference  method 
or  the  shooting  method,  and  in  the  frequency  domain  by  the  harmonic  balance 
method  [2-4]. 

The  finite-difference  method  is  not  as  widely  used  as  the  shooting  method 
since  it  requires  a  large  amount  of  memory  and  tends  to  have  more  convergence 
problems  than  the  shooting  method  [2].  The  shooting  method  is  a  popular  time- 
domain  method  in  which  a  two-point  boundary  value  problem  is  formulated.  It 
uses  the  transient  analysis  over  one  period  to  obtain  the  final  value  at  the  end  of 
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the  period,  so  it  can  be  viewed  as  an  accelerated  transient  method.  The 
acceleration  is  obtained  by  adjusting  the  initial  condition  at  the  beginning  of  each 
one-period  transient  simulation.  The  shooting  method  works  well  for  highly 
nonlinear  circuits  since  it  is  based  on  a  time-domain  analysis  and  has  the  ability  to 
handle  the  nonlinear  behavior  of  circuits.  However,  the  shooting  method  is  not 
suitable  for  circuits  with  wide  frequency  separations  and  it  cannot  handle 
distributed  devices  commonly  encountered  in  RF  simulation  [2], 

The  harmonic  balance  (HB)  method  is  a  frequency-domain  algorithm  used 
for  high  accuracy  computation  of  the  periodic  steady  state  of  RF  circuits  [2-3].  It 
is  well  suited  for  RF  system  simulation  where  the  spectral  and  dynamic  range 
requirements  necessitate  accurate  analysis  of  harmonic  content  with  large 
frequency  differences  [6].  Because  it  is  a  frequency-domain  technique,  distributed 
models  can  be  easily  and  accurately  included  in  the  simulation  [2].  With  a  Krylov 
subspace  solver,  the  harmonic  balance  method  can  be  applied  to  full-chip 
simulation  with  multi-tone  excitations  [6] . 

This  chapter  is  organized  as  follows.  The  harmonic  balance  method 
formulation  is  reviewed  in  Section  2.2.  The  discrete  Fourier  transform  for  a 
periodic  signal  is  also  presented.  Newton’s  method  for  the  harmonic  balance 
problem  is  described  in  Section  2.3.  The  matrix  implicit  Krylov  subspace  method 
is  briefly  discussed  in  Section  2.4.  The  spectral  accuracy  of  the  harmonic  balance 
method  compared  to  the  time-domain  methods  is  described  in  Section  2.5.  The 
error  mechanisms  in  the  harmonic  balance  method  are  introduced  in  Section  2.6, 
and  a  chapter  summary  is  provided  in  Section  2.7. 
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2.2  Harmonic  Balance  Equation  Formulation 

The  harmonic  balance  method  is  a  frequency-domain  method  where  the 
circuit  waveforms  are  represented  by  the  Fourier  sine  and  cosine  series  [2-3].  This 
approximation  of  a  time-domain  waveform  as  a  Fourier  series  naturally  and 
efficiently  guarantees  that  the  solution  obtained  is  indeed  the  periodic  steady  state 
of  the  system.  The  unknowns  in  the  harmonic  balance  method  are  the  frequency- 
domain  Fourier  coefficients  and  not  the  time-domain  waveforms,  whereby  the 
dynamic  aspects  of  the  problem  are  eliminated. 


2.2.1  The  Discrete  Fourier  Transform  (DFT) 

The  harmonic  balance  method  formulates  circuit  equations  in  the  frequency 
domain.  The  Fourier  transform  is  used  to  transform  the  circuit  unknowns  between 
the  time  and  frequency  domains.  Since  the  time-domain  circuit  unknowns  are 
real-valued  waveforms,  a  simplified  version  of  the  standard  discrete  Fourier 
transform  (DFT)  is  introduced  here. 

A  periodic  waveform  with  period  T  can  be  represented  by  a  Fourier  series: 

x(t)  =  X0  +  jr  X/  cos (cofj+X*  sin  (cof)  (2.1) 

i=i 

where  (d.  is  the  frequency  for  the  zth  harmonic.  Xo  is  the  DC  value  and  X ( X ) 

is  the  Fourier  coefficient  of  the  cosine  (sine)  part  of  the  zth  harmonic.  To  make  the 
problem  computationally  tractable,  it  is  necessary  to  truncate  the  frequencies  to  a 
finite  set.  If  the  energy  in  harmonics  higher  than  K  is  negligible,  the  waveform  x(t) 
can  be  approximated  as 

K 

x(t)  ~  X0  +  ^  X ■  cos(coit)+  X sin  (cof) 

i=i 


(2,2) 
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As  the  spectrum  of  x(t )  is  finite,  it  is  possible  to  sample  the  time-domain 
waveform  at  a  finite  number  of  time  points  and  calculate  the  Fourier  coefficients. 
To  get  2K+ 1  Fourier  coefficients,  at  least  the  same  number  of  sampled  time  points 
have  to  be  used.  This  results  in  a  set  of  2K+\  equations  with  2K+1  unknowns: 


where 


*0 

i 

X 

o 

_ 1 

A 

a; 
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(2.4) 


is  the  inverse  DFT,  and  x\  is  the  value  of  x(t)  sampled  at  equally  spaced  time 
points  tk( k  =  0,  1,  ...,  2 K),  where  tk  =  kT  /(2 K  +  1).  This  system  is  invertible  and 

can  be  compactly  written  as  r1  A  =  x .  Inverting  T  1  gives  Tx  -  X  .  T  (DFT) 
and  r1  (inverse  DFT)  are  a  discrete  Fourier  transform  pair. 


2.2.2  Formulating  Harmonic  Balance  Equations 

The  circuit  behavior  can  be  described  in  the  time  domain  by  a  system  of 
equations: 


10 


i(x(t))  +  —  q(x(t))  +  s(t )  =  0  (2.5) 

dt 

where  x  is  the  vector  of  circuit  waveforms  (voltage/current),  i  is  a  function  of  the 
contributions  from  the  nonreactive  elements  to  the  circuit  equations,  q 
(charge/flux)  is  the  contribution  from  reactive  elements,  and  s  is  an  input  source. 

Let  the  circuit  be  driven  by  a  single  periodic  input  source  with  period  T. 
Then  the  circuit  steady-state  response  x(t)  will  also  be  periodic  and  a  truncated 
Fourier  series  can  be  used  to  approximate  the  signal  as  in  Eq.  (2.2).  Using  the 
Fourier  coefficients  as  unknowns,  the  circuit  equations  can  be  formulated  in  the 
frequency  domain  as 

Fi(FlX)  +  QFq(  T-'X)  +  5=0  (2.6) 

where  X  is  the  vector  of  circuit  waveforms  in  the  frequency  domain  represented 
by  their  Fourier  coefficients,  F  and  F  are  the  discrete  Fourier  transform  matrix 
pair,  S  is  the  frequency-domain  representation  of  the  stimulus  vector,  Fl  is  a 
block-diagonal  matrix  representing  the  time  derivative  operation  [2] : 

0 

G7l 

where  0)k  is  the  frequency  for  the  kth  harmonic.  Since  each  frequency  requires 

two  Fourier  coefficients,  except  DC,  the  number  of  Fourier  coefficients  for  each 
circuit  waveform  is  2K+1.  If  there  are  N  waveforms  in  the  circuit,  the  total 
number  of  unknowns  is  N(2K+1)  [23-24]. 

The  harmonic  balance  equations  can  be  seen  as  KCF  formulated  in  the 
frequency  domain  where  the  frequency  spectrum  at  each  node  is  balanced.  The 
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advantage  of  this  set  of  equations  is  that  the  time  domain  derivative  operation  ( — ) 

dt 

is  changed  into  an  algebraic  multiplication  operation  (fl)  in  the  frequency 
domain.  This  directly  yields  the  solution  in  the  steady  state. 

The  harmonic  balance  method  can  also  be  written  in  the  time  domain  by 
converting  Eq.  (2.6)  from  the  frequency  domain  to  time  domain: 

i(x)  +  r~1Qrq(x)  +  s  =  0  (2.8) 

where  x  and  s  are  the  vectors  of  2K+ 1  time-domain  samples  of  unknown 
waveforms  and  stimuli,  respectively.  From  Eqs.  (2.5)  and  (2.8)  we  can  see  that  the 
time  derivative  in  Eq  (2.5)  is  approximated  by  D  =  F  <TF  in  the  harmonic 
balance  method.  The  time-domain  differentiation  matrix  D  is  a  block-dense,  real 
matrix.  The  harmonic  balance  method  can,  therefore  be  viewed  as  a  finite- 
difference  method  of  order  K ,  the  order  of  the  Fourier  truncation  [5]. 


2.3  Newton’s  Method  for  the  Harmonic  Balance  Problem 

The  harmonic  balance  method  results  in  a  system  of  nonlinear  algebraic 
equations.  Newton’s  method  is  applied  to  solve  these  equations  iteratively.  This 
application  results  in  the  iteration 

Jm(XM-X,)  =  -F(Xi)  (2.9) 

where  i  is  the  Newton  iteration  index.  JU)  is  the  ith  Newton  iteration  Jacobian 

j<l>  =rG'T  1  +nrc(,T1  (2.10) 

where  G  and  C  are  block  diagonal  matrices,  with  the  diagonal  elements 
representing  the  circuit  linearized  at  the  sampled  time  points: 
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C° 

,  c. 

C  = 


c „ 


(2.11) 


(2.12) 


where  Gk  = 


di(x ) 


'xk 


and  Ck  = 


dq(x) 

dx 


xk 


are  the  NxN  linearized  conductance  and 


capacitance  matrices  at  time  point  tk,  respectively.  These  linearized  matrices  can 
be  obtained  by  multiple  small- signal  analyses  around  the  operating  points 
specified  at  each  sample  time  point. 


For  an  arbitrary  nonlinearity,  the  Fourier  coefficients  of  the  response  can  not 
be  computed  directly  from  the  Fourier  coefficients  of  the  stimulus.  This  is 
possible  only  for  some  special  cases  [25].  Therefore,  the  nonlinear  device  is  first 
evaluated  in  the  time  domain  and  then  transformed  to  the  frequency  domain  by 
the  Fourier  transform. 


Let’s  consider  a  nonlinear  resistive  example 

i  =  i(v) 


(2.13) 


where  i  is  the  current  through  the  nonlinear  resistor  and  v  is  the  voltage  across  it. 
The  frequency  domain  equation  for  the  nonlinear  resistor  is  given  by 

/  =  ri(r-V)  (2.14) 


where  I  and  V  are  the  Fourier  coefficients  of  the  current  and  voltage,  respectively. 
The  frequency-domain  voltage  Fourier  coefficients  are  first  transformed  into  a 
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time-domain  circuit  waveform  v  =  F  lV  and  then  the  current  i(v)  is  evaluated  in 
the  time  domain.  The  frequency-domain  current  Fourier  coefficients  are  obtained 
by  the  use  of  the  Fourier  transform  V  on  the  time-domain  current  waveform,  i. 

When  using  Newton’s  method,  the  frequency-domain  current  is  linearized  at 
each  iteration 


j(M)  _  n(r~V(i+1))  =  r/(r~V(0)+r—  — _ —  \f(,+1) 

1  J  1  j  3(r"V(0)  dv(i)  ^ 

=  r(G(v(0  )v(i+I)  +  (i(v(0 )  -  G(v(0  )v(0 )) 


-y(0) 


(2.15) 


From  Eq.  (2.15)  we  can  see  that  to  evaluate  the  nonlinear  elements  in  the 
frequency  domain  via  the  time  domain,  a  time-domain  companion  network  [26] 
can  be  used.  The  only  difference  is  that  the  companion  network  needs  to  be 
transformed  to  the  frequency  domain  by  the  Fourier  transform  matrix  T . 

For  linear  circuits,  the  G  and  C  matrices  are  independent  of  the  voltages 
which  results  in  an  iteration- invariant  block-diagonal  Jacobian  matrix.  Newton’s 
method  converges  in  one  iteration,  and  the  harmonic  balance  method  is  equivalent 
to  a  phasor  analysis.  For  nonlinear  circuits,  the  Jacobian  is  no  longer  block 
diagonal.  The  off-diagonal  elements  represent  the  frequency  coupling  whose 
strength  depends  on  the  nonlinearity  of  the  circuit  [2] . 

Newton’s  method  may  only  converge  if  the  initial  guess  is  close  enough  to 
the  solution.  Thus,  a  good  initial  guess  is  important  in  determining  the  likelihood 
of  convergence.  For  many  circuits,  the  DC  operating  point  can  be  used  as  a  good 
initial  guess.  An  initial  guess  can  also  be  generated  by  linearizing  the  circuit  at 
the  DC  operating  point,  applying  the  stimulus,  and  performing  a  phasor  analysis 
at  each  frequency. 

One  feature  of  Newton’s  method  is  that,  once  a  solution  has  been  found,  the 
resulting  harmonic  Jacobian  represents  a  linearization  of  the  circuit  about  its  time 
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varying  operating  point.  This  can  be  used  in  small- signal  time  varying  analysis 
such  as  mixer  noise  analysis  or  phase  noise  calculation  for  oscillators. 


2.4  Matrix-implicit  Krylov-subspace  Methods 

The  Jacobian  matrix  of  the  harmonic  balance  method  is  a  dense  matrix  of 
size  N  (2K+X).  The  factorization  of  such  a  large  dense  matrix  by  a  direct  method 
is  prohibitively  expensive,  which  limits  the  application  of  the  harmonic  balance 
method  to  larger  circuits. 

Krylov  subspace  methods  have  been  shown  to  be  effective  when  applied  to 
the  solution  of  a  large-scale  system  of  harmonic  balance  equations.  The  Krylov 
subspace  methods  involve  only  matrix  vector  products.  Therefore  these  methods 
can  be  easily  applied  in  situations  where  the  matrix  may  not  be  directly  available, 
and  generating,  storing,  and  multiplying  with  that  matrix  involves  significant 
overhead.  For  instance,  in  the  harmonic  balance  method,  the  coefficient  matrix  is 
available  as  a  sequence  of  transforms  which  can  be  efficiently  applied  to  a  vector 
[27-31], 

The  convergence  rate  of  a  Krylov  subspace  method  strongly  depends  on  the 
condition  number  or  spectral  properties  of  the  coefficient  matrix.  The  idea  of 
preconditioning  is  to  transform  the  original  linear  system  into  one  that  is 
equivalent  with  the  same  solution,  but  has  more  favorable  spectral  properties  or 
condition  number.  The  preconditioner  is  a  matrix  that  performs  this 
transformation.  One  way  to  design  a  preconditioner  is  to  construct  a  matrix  that 
approximates  the  coefficient  matrix,  and  is  easy  to  invert.  The  successful  use  of 
Krylov  subspace  methods  in  most  situations,  is  determined  by  the  ability  to  form 
an  appropriate  preconditioner  [27-31]. 
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In  the  frequency-domain  formulation  of  the  harmonic  balance  method,  the 
most  commonly  used  preconditioner  is  a  block  diagonal  matix  which  is  formed 
from  the  Jacobian  matrix  by  only  keeping  the  linear  pail.  This  preconditioner  can 
be  easily  and  efficiently  formed  and  inverted  due  to  its  block-diagonal  structure. 
The  block-diagonal  preconditioner  has  been  shown  to  be  very  effective  for  a  large 
number  of  weakly  nonlinear  problems.  However,  it  becomes  increasingly 
ineffective  for  strongly  nonlinear  problems  because  it  becomes  a  poor 
approximation  of  the  harmonic  Jacobian  matrix  [27]. 

Adaptive  preconditioners  [28-29]  are  formed  by  neglecting  the  small 
harmonics  of  nonlinear  conductances  and  capacitances.  This  process  corresponds 
to  retaining,  in  addition  to  the  block-diagonal  entries,  only  the  most  essential  non¬ 
diagonal  entries  of  the  Jacobian  matrix. 

Recently,  the  Jacobian  matrix  from  a  lower-order  finite  difference  method 
on  a  uniform  grid  was  also  used  as  a  preconditioner.  This  was  possible  since  the 
time-domain  formulation  of  the  harmonic  balance  method  is  a  finite-difference 
method  in  disguise  [30-31].  The  Jacobian  matrix  from  a  finite-difference  method 
is  a  good  approximation  of  the  Jacobian  matrix  of  the  harmonic  balance  method 
formed  in  the  time  domain.  The  order  of  the  time-domain  integration  method 
determines  the  bands  of  the  blocks. 

With  more  complex  preconditioners,  one  can  expect  the  computational 
effort  to  reduce  due  to  a  decrease  in  the  number  of  linear  iterations.  However, 
there  will  be  an  increased  effort  in  factoring  the  preconditioners. 


2.5  Spectral  Accuracy  of  the  Harmonic  Balance  Method 

Time-domain  methods  (finite-difference  method  and  shooting  method) 
compute  the  steady  state  by  first  discretizing  the  solution  domain  [0,  7].  The  time 
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derivatives  arc  approximated  with  finite-order  integration  formulas.  For  example, 
for  the  A'-th-order  BDF  formula,  the  time  derivative  of  a  function  is  computed 
using  values  of  the  function  at  K  nearby  time  points  [2],  For  this  reason,  the 
approximation  of  the  derivative  in  the  time  domain  captures  only  the  local 
behavior  of  the  function.  These  time-domain  methods  are  exact  for  polynomials 
of  order  <  K.  In  other  words,  the  solution  is  represented  by  a  sequence  of  low- 
order  polynomials  connected  at  the  discrete  time-points.  These  methods  can, 

is 

therefore,  at  best  achieve  polynomial  convergence  with  global  errors,  i.e.,  0(h  ) 

[5]. 

The  frequency-domain  harmonic  balance  method  approximates  the  solution 
waveforms  using  a  weighed  Fourier  basis.  The  time-derivative  of  this 
approximation  can  be  computed  exactly  as 

K 

x(t)  =  'y't(a)iXis  cos (cOft)- coiX ^  sin (coj))  (2.16) 

i=l 

If  the  circuit  waveform  is  sampled  at  2 K+l  uniformly  distributed  time  points, 
then  the  time  derivative  at  each  time  point  is  approximated  with 

*<>  i  r  *o 

xx  xx 

x 2  =r~1ar  x?  (2.17) 

X2K- 1  X2K- 1 

X2K  J  L  X2k 

Since  r_1£2r  is  a  dense  matrix,  the  time  derivative  of  the  circuit  waveform  at 
each  time  point  is  computed  using  values  of  the  circuit  waveforms  at  2 K+l  time 
points.  Since  these  time  points  span  the  whole  period  of  the  circuit  waveform,  the 
time  derivative  in  the  harmonic  balance  method  includes  the  global  behavior  of  a 
function. 
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Due  to  the  global  nature  of  the  harmonic  balance  method,  as  i  — >  °o  the  zth 
Fourier  coefficient  X,  decays  rapidly  if  the  solution  waveform  is  smooth.  This 
rapid  decay  of  the  coefficients  means  that  the  truncated  Fourier  approximation  of 
a  solution  expanded  using  a  few  additional  terms  (by  slightly  increasing  K) 
represents  an  exceedingly  good  approximation  of  the  solution,  a  property  known 
as  spectral  accuracy,  exponential  convergence,  or  infinite-order  accuracy.  The 
harmonic  balance  method  will  not  achieve  spectral  accuracy  if  a  discontinuity  is 
present  in  the  solution  waveforms  or  one  of  their  derivatives.  If  the  device 
models  are  p- times  continuously  differentiable  functions  and  periodic  in  all  its 
derivatives,  then  the  errors  for  the  harmonic  balance  method  will  be  0(KP )  [5]. 


2.6  Error  Mechanisms 


The  errors  in  the  harmonic  balance  method  come  from  the  truncation  of  the 
frequency  spectrum  to  some  finite  order  and  the  use  of  a  finite-order  discrete 
Fourier  transform.  This  error  can  be  split  into  two  parts,  a  truncation  error  and  an 
aliasing  error  [2]. 

Consider  a  circuit  with  periodic  or  quasi-periodic  solutions.  Let  F  be  the 
sum  of  currents  at  every  node  and  every  frequency.  is  defined  as  the  truncated 
set  of  frequencies.  Let  Fexact  be  the  result  when  the  circuit  equations  are  evaluated 
without  an  error  at  all  frequencies.  Then  the  truncation  error  is  defined  as 


K»jv,k) 


fO  if  cok  e  Ak 

1  Fe*,«(v’k)  otherwise. 


(2.18) 


Thus,  Ftnmc  represents  all  currents  generated  by  the  circuit  at  frequencies  other 
than  those  in  [2]. 
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Aliasing  errors  are  introduced  when  a  time-domain  waveform  is  converted 
into  a  frequency-domain  spectrum  using  a  finite-order  discrete  Fourier  transform. 
Aliasing  is  caused  by  using  too  few  time  points  to  represent  the  time  waveform. 
The  Nyquist  criterion  states  that  to  avoid  aliasing  errors,  the  sampling  frequency 
should  be  at  least  twice  the  highest  frequency  present  in  the  time-domain 
waveform  [32].  If  not,  an  error  is  introduced  in  the  low-frequency  Fourier 
coefficients.  In  this  sense  the  harmonic  balance  method  is  more  accurate  for  band- 
limited  signals  where  the  error  from  aliasing  does  not  corrupt  the  frequency 
spectrum.  The  only  way  to  eliminate  aliasing  errors  is  to  increase  the  number  of 
data  points  by  oversampling  the  time-domain  waveform. 

The  aliasing  error  is  generally  as  large  as  or  larger  than  the  truncation  error 
[2],  One  approach  to  estimate  the  error  due  to  the  truncation  and  aliasing  is  to 
simply  increase  the  number  of  frequencies  used  in  the  Fourier  series 
representation  of  the  circuit  unknowns  and  re- simulate  the  circuit.  However,  this 
method  of  estimating  the  error  is  an  inexpensive  way  of  determining  whether 
enough  frequencies  have  been  included  in  the  harmonic  balance  analysis  after  a 
solution  has  been  computed. 


2.7  Summary 

This  chapter  first  presented  an  overview  of  the  harmonic  balance  algorithm. 
The  harmonic  balance  method  overcomes  several  difficulties  faced  by  time- 
domain  methods  by  solving  the  system  of  equations  in  the  frequency  domain.  It 
also  provides  higher  spectral  accuracy  compared  to  the  time-domain  methods. 
The  formulation  of  the  harmonic  balance  equations  and  solution  with  Newton’s 
method  is  described.  Error  mechanisms  in  the  harmonic  balance  method  have  also 
been  discussed. 
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3.  SIMULATION  OF  OSCILLATORS  WITH  THE  HARMONIC 

BALANCE  METHOD 


3.1  Introduction 

Oscillators  are  widely  used  in  many  applications  such  as  phase-locked  loops 
(PLL)  and  frequency  synthesizers  to  generate  reference  signals  [1],  Phase  noise  is 
an  important  characteristic  of  oscillators  and  can  be  determined  only  if  the  circuit 
operates  in  its  periodic  steady  state  [33-34].  For  this  reason,  periodic  steady-state 
simulations  of  oscillators  are  extremely  important  [2-3]. 

The  periodic  steady- state  analysis  can  be  done  in  either  the  time  domain  or 
in  the  frequency  domain.  The  frequency-domain  harmonic  balance  analysis  is 
spectrally  more  accurate  than  the  time-domain  method  since  low-order  integration 
is  typically  used  in  the  latter.  It  is  also  very  efficient  for  RF  system  simulation 
since  it  can  handle  the  wide  spread  of  frequencies  naturally. 

The  harmonic  balance  method  has  been  applied  to  oscillator  simulation  by 
adding  the  unknown  oscillation  frequency  to  the  set  of  circuit  state  variables. 
However,  oscillator  simulation  using  the  harmonic  balance  method  has  proven  to 
be  difficult  due  to  a  small  region  of  convergence  and  the  existence  of  the 
degenerate  DC  solution  [2].  Careful  implementation  and  special  techniques  must 
be  used  for  convergence  [2,  7]. 

In  this  chapter,  the  basic  oscillator  structures  and  their  operation  are  first 
introduced  in  Section  3.2.  The  application  of  the  harmonic  balance  method  to 
oscillator  simulation  is  reviewed  in  Section  3.3.  The  probe  technique  for 
improving  the  convergence  in  frequency-domain  oscillator  simulation  is  described 
in  Section  3.4. 
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3.2  Oscillators 

Most  electronic  signal  processing  systems  require  frequency  or  time 
reference  signals.  In  communication  systems  oscillators  generate  the  signals  to 
modulate  the  baseband  message  signal  in  transmitters,  and  process  the  received 
signals  in  receivers.  There  are  many  different  kinds  of  oscillators  for  various 
applications. 


3.2.1  Oscillator  Structures 

Oscillators  can  be  broadly  classified  into  three  categories:  harmonic 
oscillators,  ring  oscillators,  and  relaxation  oscillators  [1]. 

The  most  common  harmonic  oscillators  are  those  that  use  resonant  circuits 
consisting  of  inductors  and  capacitors.  Harmonic  oscillators  that  use  LC  resonant 
circuits  are  known  as  LC  oscillators.  A  typical  LC  oscillator  is  made  up  of  three 
distinct  blocks  as  shown  in  Fig.  3.1.  The  amplifier  supplies  energy  to  maintain 
oscillations  in  the  circuit.  The  LC  tank  circuit  alternately  stores  energy  in  the 
inductor  and  the  capacitor.  A  portion  of  the  output  of  the  LC  network  is  fed  back 
to  the  input  of  the  amplifier  through  a  positive  feedback  network.  This  helps  to 
sustain  oscillations  by  overcoming  the  effect  of  the  damping  caused  by  the 
internal  and  load  resistances  of  LC  tanks. 
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Figure  3.1:  Block  diagram  of  an  LC  oscillator. 


The  operating  or  center  frequency  (/)  of  the  LC  oscillator  is  fixed  by  the 
resonant  frequency  of  the  LC  tank  circuit 


j  ~ - 1= 

27r4LC 

An  important  quantity  in  the  characterization  of  LC  tank  circuits  is  the 
quality  factor,  Q ,  which  is  defined  as 


Q  =  a) 


energy  stored 
average  power  dissipated 


The  resonator  Q  strongly  influences  both  the  phase  noise  and  power  consumption 
of  an  oscillator.  The  inductor  in  an  LC  oscillator  is  usually  the  most  critical  circuit 
element  in  the  design.  Typically,  the  Q  of  the  inductor  dominates  the  total  Q  of 
the  tank  circuit. 

Ring  oscillators  are  resonator-less  oscillators  consisting  of  an  odd  number  of 
single-ended  inverters  or  an  even/odd  number  of  differential  inverters  with 
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appropriate  connections.  A  single-ended  ring  oscillator  structure  is  shown  in  Fig. 
3.2.  If  each  inverter  is  ideal  this  circuit  is  a  form  of  negative  feedback,  but  in 
reality  each  inverter  exhibits  a  lagging  phase  shift  arising  from  the  propagation 
delay.  Enough  lagging  phase  shift  (at  least  three  inverters  in  the  loop)  added  to  the 
inversions  ultimately  turns  the  negative  feedback  into  a  positive  feedback.  As  a 
result  the  circuit  is  unstable  and  oscillations  occur  [35].  In  each  half-period,  the 
signal  propagates  around  the  loop  with  an  inversion.  Whereby  the  oscillation 
frequency  for  a  ring  oscillator  is  determined  by 


/  = 


1 

2  nrd 


(3.3) 


where  n  is  the  number  of  stages  in  the  ring  oscillator,  and  Td  is  the  time  delay  of 
one  stage  [35]. 


Figure  3.2:  Block  diagram  of  a  ring  oscillator. 


The  operation  of  a  relaxation  oscillator  is  very  different  from  the  other  types 
of  oscillators.  Relaxation  oscillators  have  poles  located  in  the  right  half  plane  to 
make  it  an  unstable  circuit.  Unlike  LC  oscillators,  the  poles  for  relaxation 
oscillators  are  real  instead  of  complex  as  shown  in  Figure  3.3.  Therefore,  the 
signal  amplitude  will  grow  exponentially  until  some  transistors  are  cut  off  [36]. 
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(a)  (b) 

Figure  3.3:  (a)  Pole  location  for  relaxation  oscillators,  (b)  Exponential  growth  of 
circuit  waveform  for  a  relaxation  oscillator. 

The  relaxation  oscillator  has  two  distinct  phases  of  operation,  a  fast 
regenerative  switching  phase  and  a  recovery  phase.  Oscillations  in  relaxation 
oscillators  are  commonly  achieved  by  the  charging  and  discharging  of  capacitors 
with  a  constant  current  between  two  threshold  levels.  The  oscillation  frequency  of 
a  relaxation  oscillator  is  set  by  the  charge/discharge  rate  of  the  capacitor  and  the 
difference  between  two  separate  threshold  levels  at  the  input  of  the  regenerative 
element. 

Oscillators  employed  in  RF  applications  must  meet  stringent  phase  noise 
requirements,  so  LC  type  oscillators  with  a  high  quality  factor  are  preferred. 
High-Q  factors  cannot  be  achieved  in  integrated  implementations  because  of  the 
availability  of  low  Q  inductors.  The  trend  toward  large- scale  integration  and  low 
cost  makes  it  desirable  to  implement  oscillators  on  chip.  Ring  and  relaxation 
oscillators  have  raised  interest  among  circuit  designers  because  they  are  very 
compact,  simple  and  can  operate  at  high  speeds  [1].  A  ring  oscillator  can  also 
readily  yield  output  phases  in  quadrature. 
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3.2.2  Oscillator  Analysis 

To  generate  a  periodic  output,  the  oscillator  circuit  must  have  a  self- 
sustaining  mechanism  that  allows  its  own  noise  to  grow  and  eventually  become  a 
periodic  signal.  There  are  two  analysis  points  of  view  for  oscillators,  the  negative 
resistance  approach  and  the  feedback  approach  [36]. 

The  negative  resistance  concept  is  illustrated  in  Fig.  3.4.  A  RLC  tank  circuit 
has  losses  represented  by  the  resistance  R ,  so  it  cannot  oscillate  indefinitely  as  the 
stored  energy  is  dissipated  in  the  resistor  for  every  cycle.  To  sustain  oscillations 
the  loss  of  energy  must  be  compensated.  An  active  network  is  used  to  generate  a 
resistance  equal  to  - R .  In  this  case,  the  circuit  can  be  seen  as  a  lossless  LC  tank 
and  the  energy  lost  in  R  is  replenished  by  the  active  circuit  in  every  cycle. 


Rin—R 


L 


Active 

Circuit 


Figure  3.4:  Negative  resistance  approach  for  oscillator  analysis. 


The  feedback  approach  is  illustrated  in  Fig.  3.5.  The  oscillator  is  viewed  as 
a  circuit  with  positive  feedback  where  a  portion  of  the  output  is  combined  in 
phase  with  the  input.  The  closed-loop  gain  is  given  by 

X0  _  Ajs) 


Xs  1  -A(s)B(s) 


(3.4) 
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where  A(s)B(s )  is  the  loop  gain.  The  oscillation  reaches  its  steady  state  when  5 
is  purely  imaginary  and 

ZA(ja>Q)  + ZB(jo)0)  =  2mi  where/?  =  1, 2, 3- 
where  (Oo  is  the  oscillation  frequency  [1]. 


3.3  Application  of  the  Harmonic  Balance  Method  for  Oscillator  Simulation 

An  important  application  of  the  harmonic  balance  method  is  in  determining 
the  steady-state  behavior  of  oscillators  [2].  Oscillators  present  two  problems  not 
found  in  forced  circuits.  First,  the  period  of  oscillation  is  an  unknown  and  must  be 
determined.  Second,  there  is  no  input  to  fix  the  phase.  Thus,  if  one  solution  exists, 
then  an  infinite  continuum  of  phase-shifted  solutions  exists.  For  these  reasons,  Eq. 
(2.6)  is  modified  to  handle  oscillators  by  adding  the  fundamental  frequency  to  the 
list  of  unknowns.  An  additional  equation  is  also  added  to  fix  the  phase.  The 
common  practice  is  to  set  the  sine  part  of  the  fundamental  of  the  chosen  signal  to 
zero  [2]. 

F(x,co)  =  n(rlx)  +  Q(co)rq(rlx)  +  s  =  o 

*:(D= 0 


(3.6) 
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This  set  of  equations  can  be  solved  using  Newton’s  method  where  the 
Jacobian  is  given  by 


J 


osc 


J 

«(1))7 


a f 

dco 

0 


(3.7) 


with 


dF  _  nrqcr-'x) 

dco  CO 


(3.8) 


where  J  is  the  Jacobian  in  Eq.  (2.10)  and  e'm  (1)  is  the  unit  vector  that  selects  the 
sine  part  of  the  first  harmonic  of  the  chosen  signal. 

This  implementation  is  straightforward.  However,  without  an  extremely 
good  initial  guess  for  both  the  oscillation  frequency  and  all  the  circuit  unknowns, 
the  method  either  tends  to  converge  to  the  trivial  DC  solution,  or  fails  to  converge 
[7-9].  The  reason  for  the  convergence  difficulty  is  that  the  unknown  frequency  is 
adjusted  at  every  iteration  along  with  the  node  voltages.  The  node  voltages  at  the 
intermediate  iterations  do  not  satisfy  KCL.  Updating  the  frequency  based  on  these 
node  voltages  can  lead  to  divergence,  or  make  convergence  difficult. 


3.4  The  Probe  Technique 

To  handle  the  convergence  problems  encountered  in  frequency-domain 
oscillator  simulation,  Ngoya  et  al.  [7]  presented  the  idea  of  a  probe  which  has 
become  a  standard  technique  for  frequency-domain  oscillator  simulation.  The 
probe  can  be  a  voltage  probe  or  a  current  probe  as  shown  in  Fig.  3.6.  A  probe  is 
composed  of  an  ideal  independent  source  (current  source  or  voltage  source  for 
current  or  voltage  probe,  respectively)  associated  with  a  filter.  The  filter 
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electrically  disconnects  the  probe  from  its  termination  for  any  frequency  different 
from  that  of  the  probe  source  (a  short  circuit  for  the  current  probe,  and  an  open 
circuit  for  the  voltage  probe). 


(a) 


(b) 


Figure  3.6:  (a)  Voltage  probe,  (b)  Current  probe. 
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To  simulate  an  oscillator  the  voltage  probe  is  connected  in  parallel  and  the 
current  probe  is  connected  in  series.  The  implementation  of  the  current  probe  is 
not  as  straightforward  as  that  of  the  voltage  probe  since  the  circuit  must  be  broken 
at  a  certain  point  to  insert  the  current  probe.  Furthermore,  it’s  easier  to  specify  an 
initial  guess  for  a  voltage  than  a  current,  so  the  voltage  probe  is  more  commonly 
used. 

The  simulation  of  an  oscillator  with  the  aid  of  the  voltage  probe  is  shown  in 
Fig.  3.7  where  the  probe  voltage  and  frequency  keep  changing  until  the  solution  is 
obtained  when  the  probe  current  becomes  zero.  Then  all  the  circuit  equations  are 
satisfied  without  the  probe  being  a  part  of  the  circuit.  In  this  case,  the  voltage 
probe  is  actually  detached  from  the  oscillator  circuit.  The  advantage  of  the  probe 
approach  is  that  the  simulation  of  an  oscillator  circuit  is  changed  into  a  set  of 
closely  related  forced  circuits,  which  can  be  easily  handled  using  the  standard 
harmonic  balance  method. 


Figure  3.7:  Oscillator  simulation  with  the  aid  of  the  voltage  probe. 
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The  node  at  which  the  voltage  probe  is  attached  in  the  oscillator  circuit 
influences  the  robustness  of  the  approach.  It  is  observed  that  better  convergence  is 
achieved  if  the  voltage  probe  is  attached  at  a  high  impendence  node.  The  collector 
node  of  the  BJT  device  or  the  drain  node  of  the  MOS  device  is  usually  a  good 
choice. 


3.4.1  The  Two-level  Method 


The  voltage  probe  is  implemented  by  a  two-level  Newton  method  [7].  The 
bottom  level  of  the  two-level  method  solves  the  harmonic  balance  equations  of 
the  voltage-probe  forced  circuit.  The  voltage  probe  is  implemented  directly  in  the 
frequency  domain  with  the  corresponding  branch  constitutive  equations.  The  sine 
part  of  the  probe  voltage  is  set  to  zero  to  fix  the  phase  and  thus  selects  one  of  the 
many  equivalent  solutions.  The  linear  equations  solved  at  each  step  of  the  bottom- 
level  Newton  method  are  given  by 
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(3.9) 


where  J  is  the  Jacobian  in  Eq.  (2.10).  ecm(l)  (V  (1))  is  the  unit  vector  that  selects 
the  cosine  (sine)  part  of  the  first  harmonic  of  the  chosen  node  to  which  the 
probe  is  attached ■  Icprobe(Isprobe)  is  the  cosine  (sine)  part  of  the  probe  current,  RhsF 

is  the  right-hand-side  vector  of  the  circuit  which  includes  the  circuit  stimulus  and 
the  contributions  from  nonlinear  devices. 

After  getting  the  probe  current  from  the  bottom  solution,  the  top-level  solves 
the  probe  equations  as  shown  in  Eq.  (3.10)  to  update  the  probe  voltage  and  the 
frequency. 
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ICprobXVprobe,0))=0 

lloJVPro^  ®)  =  0 


(3.10) 


The  updated  frequency  and  probe  voltage  are  then  applied  to  solve  the  bottom- 
level  probe  forced  circuit  again.  The  process  is  repeated  until  the  probe  current 
goes  to  zero  where  the  solution  is  reached. 


The  Jacobian  for  the  top-level  probe  equations  is: 
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(3.11) 


It  should  be  noted  that  the  frequency  is  updated  only  at  the  top  level  where 
the  bottom-level  probe-forced  circuit  already  converges  and  KCLs  are  satisfied 
for  the  circuit-probe  combination.  This  results  in  an  improved  convergence. 


3.4.2  Convergence  Conditions  and  Challenges 

The  success  of  the  two-level  method,  however,  depends  on  1)  a  close  initial 
guess  for  the  oscillation  frequency  and  the  probe  voltage,  and  2)  the  convergence 
of  the  bottom- level  probe-forced  circuit  equations. 

A  pole-zero  analysis  can  be  used  for  the  initial  frequency  guess.  A  curve 
tracking  is  used  to  get  the  initial  guess  for  the  probe  voltage  as  shown  in  Fig.  3.8 
where  the  curve  of  the  probe  current  is  tracked  as  a  function  of  the  probe  voltage 
for  a  fixed  frequency  [7].  If  the  frequency  used  is  the  actual  oscillation  frequency, 
then  at  a  certain  probe  voltage,  the  probe  current  will  go  to  zero  as  expected.  This 
voltage  corresponds  to  the  solution  for  the  probe  voltage.  If  the  frequency  used  is 
close  to  the  actual  oscillation  frequency,  such  as  the  frequency  from  a  pole-zero 
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analysis,  then  there  will  be  a  local  minimum.  The  corresponding  voltage  value  is 
quite  close  to  the  probe  voltage  solution  as  shown  in  Fig.  3.8  and  can  be  used  as  a 
good  initial  guess  for  the  probe  voltage.  However,  this  scheme  for  finding  an 
initial  guess  for  the  probe  voltage  fails  for  some  high-Q  oscillators  as  will  be 
shown  later. 


Figure  3.8:  Probe  current  tracking  curves  for  selecting  an  initial  probe  voltage. 
Use  of  the  frequency  from  a  pole-zero  analysis  provides  a  good  initial  guess  for 
the  probe  voltage  (minimum  in  the  tracking  curve). 


With  a  good  initial  guess  for  the  frequency  and  probe  voltage,  Newton’s 
method  can  be  applied  to  solve  the  top-level  probe  equations.  However,  the 
construction  of  the  top-level  Jacobian  matrix  needs  the  information  from  the 
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bottom-level  Jacobian  when  the  bottom-level  equations  have  converged. 
Convergence  of  the  bottom-level  voltage-probe  forced  circuit  can  be  a  challenge 
for  some  highly  nonlinear  oscillators  such  as  ring  oscillators  as  shown  later. 


3.5  Summary 

In  this  chapter  various  oscillator  structures  and  their  operation  are  reviewed. 
Harmonic  oscillators,  ring  oscillators  and  relaxation  oscillators  are  discussed.  The 
application  of  the  frequency-domain  harmonic  balance  method  to  oscillator 
simulation  is  introduced.  The  direct  method  of  simply  adding  the  frequency  to  the 
list  of  unknowns  suffers  from  convergence  problems.  To  overcome  the 
convergence  problems,  a  more  robust  probe  technique  is  typically  used.  This 
technique  is  described  along  with  its  implementation  by  a  two-level  method. 
Conditions  for  convergence  are  also  discussed  in  detail. 
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4.  FREQUENCY-DOMAIN  SIMULATION  OF  HIGH-Q  OSCILLATORS 

4.1  Introduction 

Traditional  time-domain  steady  state  analysis  is  impractical  for  high-Q 
oscillators  since  a  very  large  number  of  oscillator  periods  have  to  be  simulated 
before  the  steady  state  is  reached.  Accelerated  steady  state  analysis  methods  (in 
time  or  frequency  domains)  alleviate  this  problem  [2]. 

Frequency-domain  analysis  using  the  harmonic  balance  method  is 
particularly  efficient  for  high-Q  oscillators.  The  steady-state  circuit  waveforms 
can  be  directly  obtained  using  a  small  number  of  harmonics  because  of  the  high-Q 
nature  of  the  circuit.  Furthermore,  the  harmonic  balance  method  is  more  accurate 
than  the  time  domain  method  for  circuits  with  band-limited  waveforms  due  to  its 
exponential  convergence  nature.  Efficiency  and  accuracy  requirements  make  the 
frequency-domain  harmonic  balance  method  an  ideal  candidate  for  high-Q 
oscillator  simulation. 

The  harmonic  balance  method  has  been  successfully  used  for  the  simulation 
of  oscillators  using  the  concept  of  a  periodic  voltage  probe  [7].  However,  high-Q 
oscillator  simulation  using  the  harmonic  balance  method  is  a  difficult  problem 
due  to  a  small  region  of  convergence.  Prior  research  has  focused  on  obtaining  the 
region  of  convergence  with  two  dimensional  optimization  procedures  [37]. 
However,  this  procedure  can  be  very  expensive  since  it  involves  repeatedly 
solving  many  harmonic  balance  equations. 

A  detailed  analysis  in  [8,  9]  shows  that  a  continuation  method  in  conjunction 
with  the  voltage  probe  of  [7]  can  be  successfully  used  to  simulate  a  high-Q 
oscillator.  A  multistage  continuation  approach  for  oscillator  simulation  using  a 
Hopf  bifurcation  is  also  proposed  in  [10].  However,  the  specific  details  of  the 
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continuation  methods  are  not  provided  [8-10].  In  this  chapter,  we  explore  the  use 
of  homotopy  methods  [11-12]  for  the  steady-state  simulation  of  high-Q  oscillators 
in  the  frequency  domain. 

This  chapter  is  organized  as  follows.  The  challenges  for  high-Q  oscillator 
simulation  in  the  frequency  domain  are  introduced  in  Section  4.2.  In  Section  4.3 
the  theory  of  homotopy  methods  is  reviewed.  Then  the  homotopy  method  solver 
HOMPACK  [11]  is  described  in  Section  4.4.  Homotopy  methods  are  applied  to 
the  frequency-domain  simulation  of  high-Q  oscillators  and  different  embedding 
methods  are  given  in  Section  4.5.  This  chapter  is  summarized  in  Section  4.6. 


4.2  Convergence  Problems  in  High-Q  Oscillators 


As  shown  in  the  previous  chapter,  the  two-level  method  requires  a  good 
starting  point  for  the  oscillation  frequency  and  the  probe  voltage.  A  pole-zero 
analysis  is  commonly  used  for  estimating  the  oscillation  frequency.  With  this 
initial  frequency,  the  probe  current  is  tracked  while  stepping  the  probe  voltage. 
For  most  circuits  a  local  minimum  is  obtained  and  the  corresponding  probe 
voltage  is  then  a  good  initial  guess  [7].  However,  for  high-Q  oscillators  such  a 
minimum  does  not  exist  even  when  the  initial  oscillation  frequency  from  a  pole- 
zero  analysis  is  sufficiently  close  to  the  exact  solution  as  shown  in  Fig.  4.1.  To  see 
how  small  the  region  of  convergence  for  a  high-Q  oscillator  is,  the  probe  current 
is  tracked  for  different  frequencies  which  are  close  to  the  actual  oscillation 
frequency.  From  Fig.  4.1  we  can  see  that  a  local  minimum  of  the  probe  current  is 


obtained  only  for  frequency  initial  guesses /that  satisfy 


Ho 

/o 


<  0.005% ,  where 


/0  is  the  oscillation  frequency.  This  indicates  that  the  region  of  convergence  is 
extremely  small  and  it  is  almost  impossible  for  a  user  to  specify  a  proper  initial 
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guess  for  the  frequency  of  oscillation.  Without  a  proper  initial  guess  for  the  probe 
frequency,  an  initial  guess  for  the  probe  voltage  cannot  be  determined  and  the 
two-level  method  will  have  convergence  problems.  For  this  reason,  a  method 
which  does  not  require  an  initial  guess  for  the  probe  voltage  is  needed  for  high-Q 
oscillator  circuits.  Homotopy  methods  with  a  proper  embedding  belong  to  this 
category  as  described  in  the  next  section. 


Figure  4.1:  Probe  current  tracking  curves,  fo  is  actual  oscillation  frequency ,fpz  is 
the  frequency  from  the  pole-zero  analysis,  and  /  is  a  guess  for  the  oscillation 
frequency. 


4.3  Theory  of  Homotopy  Methods 

Homotopy  methods  belong  to  the  class  of  continuation  methods  where  the 
homotopy  parameter  does  not  necessarily  vary  monotonically  as  the  path  from  the 
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original  solution  to  the  final  solution  is  followed.  This  is  accomplished  by  making 
the  continuation  parameter  a  function  of  the  path  length  [38].  The  path  length 
always  increases  while  the  continuation  parameter  may  vary  arbitrarily. 

Consider  the  solution  of  a  system  of  n  nonlinear  equations  in  n  variables 

F(x)  =  0  (4.1) 

where  F:  R"—>R"  is  a  smooth  mapping.  A  homotopy  function  H  (x.A)  is  created 
by  embedding  a  parameter  X  which  results  in  a  system  of  equations  of  higher 
dimension 

H(xJ l)  =  0  (4.2) 

where  H:  R"  x  R—>R".  For  A  =  0,  H(x,  0)  =  G(x)  =  0  is  an  easy  problem  to  solve, 
whereas  for  2  =  1,  H(x,  1)  =  F(x)  =  0  yields  the  solution  of  the  original  problem. 
The  parameter  X  is  called  the  homotopy  parameter.  By  following  the  solutions  of 
H(x,A)  =  0  as  X  varies  from  0  to  1,  the  solution  to  F(x)  =  0  is  found  and  the 
solution  trace  is  known  as  the  zero  curve. 

The  benefit  of  this  approach  is  that  the  solution  at  a  previous  A  point  can  be 
used  as  a  good  initial  guess  for  the  current  A  point.  With  this  good  initial  guess 
only  a  few  steps  of  Newton’s  method  are  needed  to  find  the  solutions  for  each  A. 
Provided  that  H(x,  0)  =  0  has  a  unique  solution,  the  solution  x  is  bounded  along 
the  zero  curve  and  the  homotopy  path  is  bifurcation  free,  the  homotopy  method 
guarantees  a  reliable  solution  from  an  arbitrary  starting  point  [11].  This 
characteristic  makes  homotopy  methods  extremely  useful  when  an  initial  guess, 
required  for  a  local  Newton  method,  is  not  available  or  simply  is  too  difficult  to 
obtain. 

The  choice  of  the  embedding  function  H  (x.A)  is  critical  for  the  homotopy 
method  to  work  robustly  and  efficiently  [11-12],  There  are  two  general  ways  to 
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embed  a  homotopy  parameter.  The  first  method  is  a  natural  parameter  homotopy, 
where  a  system  parameter  is  chosen  and  swept  across  a  range  of  values,  and, 
therefore,  plays  the  role  of  the  homotopy  parameter  X.  The  most  commonly  used 
natural  parameter  homotopy  in  circuit  simulation  is  the  “source  stepping”. 
However,  natural  homotopies  are  often  prone  to  bifurcations  of  the  homotopy 
path.  The  second  method  is  an  artificial  parameter  homotopy,  where  the 
homotopy  parameter  X  is  embedded  in  the  nonlinear  equations  and  is  not  related 
to  any  system  parameters.  By  proper  construction  of  the  artificial  parameter 
homotopies,  bifurcations  and  other  singular  and  ill-conditioned  behavior  can  be 
avoided. 


4.3.1  Homotopy  methods  for  DC  operating  point  problem 

Homotopy  methods  have  been  used  extensively  for  finding  the  dc  operating 
points  of  “difficult-to- solve”  nonlinear  circuits  [12-15].  In  these  applications, 
variations  of  homotopy  methods  have  been  constructed  for  the  nonlinear 
equations  describing  transistor  circuits  or  device  models. 

■  Fixed-point  homotopy 

The  fixed-point  homotopy  is  based  on  the  equation 

H{x,X)  =  AF{x)-(\-X)G(x-a)  (4.3) 

where  a  new  parameter  Ge  R"  xR"  is  embedded.  This  homotopy  represents  an 
augmented  circuit  derived  from  the  original  circuit.  A  branch  consisting  of  a 

conductance  - — —G  in  series  with  a  voltage  source  au  is  connected  to  every  node. 


At  A  =  0  the  conductance  goes  to  infinity  and  the  voltage  at  each  node  is  simply 
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forced  by  a*.  When  X  =  1 ,  the  conductance  equals  to  zero,  the  added  branches  are 
disconnected  from  the  circuit  and  the  original  circuit  is  recovered. 

■  Variable-stimulus  homotopy 

The  variable- stimulus  homotopy  is  based  on  the  equation 

H(x,  X)  =  F(x,  X)-(l-  X )G(x  -  a)  (4.4) 

In  this  homotopy  mapping,  the  node  voltages  of  the  nonlinear  elements  are 
multiplied  by  X.  So  the  homotopy  starts  from  a  linear  circuit  by  setting  all  the 
voltages  across  the  nonlinear  elements  to  zero.  The  advantage  of  this  embedding 
is  that  numerical  overflows  often  encountered  in  exponential  nonlinearities  in  p-n 
junctions  are  avoided.  This  homotopy  converges  faster  to  the  desired  solution  than 
the  fixed-point  homotopy  [12-13]. 


■  Variable-gain  homotopy 

The  variable-gain  homotopy  is  designed  specially  for  bipolar  circuits,  and  is 
based  on  the  equation 

H(x,  X)  =  F( x,  X a)  -  (1  -  X )G(x  -  a)  (4.5) 

where  a  is  a  column  vector  consisting  of  transistor  forward  and  reverse  current 
gains.  When  X  =  0  the  transistor  current  gain  goes  to  zero,  so  the  circuit  only 
consists  of  resistors  and  diodes.  Such  a  circuit  has  a  unique  and  easy-to-solve  dc 
operating  point.  The  variable  gain  homotopy  converges  very  fast  for  bipolar 
circuits  [12-13]. 


The  BLHOM  homotopy 
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The  BLHOM  homotopy  is  designed  specially  for  MOSFET  circuits  [14].  A 
key  feature  of  this  homotopy  is  that  it  uses  two  X  parameters  which  are  directly 
embedded  into  the  MOSFET  level  one  model 

Lis  =  P\ygs  ivgb ,  vdb ,  vsb ,  A2  )fh(vdb  -  vsb ,  Ax )  (4.6) 

Ax  influences  the  drain-source  driving-point  characteristic  whereas  A -2  controls 
the  transfer  characteristics.  At  (Al,A2)  =  (0, 0),  there  is  no  transfer  characteristic 
(varying  vss  does  not  alter  Ids),  and  the  driving  point  characteristic  is  less  sharper 

than  that  for  the  original  MOSFET.  The  homotopy  mapping  starts  from 
(Al,A2)  =  (0, 0)  for  which  each  MOS  device  becomes  a  two-terminal  almost 
linear  resistor,  hence  the  circuit  becomes  easy  to  solve.  This  embedding  method 
works  well  for  large-scale  MOSFET  circuits. 


4.3.2  Homotopy  methods  for  time-domain  steady-state  analysis  of  oscillators 

Recently,  homotopy  methods  have  also  been  applied  to  the  steady-state 
analysis  of  oscillators  in  the  time  domain  combined  with  the  shooting  method  [16- 
17]  and  the  finite-difference  method  [18]. 

■  Homotopy  with  shooting  method 

The  original  problem  is 

J  x(0)-4>(x(0),0,7’))  =  0 
|  x,  (0)  —a  =  0 

where  x  is  the  circuit  state  vector,  including  the  voltages  across  the  capacitors  and 
currents  through  the  inductors.  <E>  is  the  state  transition  function  which  can  be 
obtained  by  a  regular  transient  analysis  starting  from  the  initial  state  v(0). 
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xj  (0)  -  a  =  0  is  an  additional  constraint  for  an  oscillator  circuit  to  isolate  the 
solution. 

The  homotopy  map  is  based  on  the  equation 

|  Mx(0)  -  <P(x(0),  0,  T))  +  (1  -  M)(x(Q)  -  xinil )  =  0  8) 

|  Mx,  (0)  -  o)  +  (1  -  X)(T  -  Tinit )  =  0 

where  xinit  is  the  initial  value  for  the  state  variables  and  Tinir  is  the  initial  period 

value  from  a  pole-zero  analysis.  By  sweeping  X  from  0  to  1,  the  problem  goes 
from  an  easy  to  solve  initial  value  problem  to  the  original  periodic  boundary  value 
problem. 


■  Homotopy  with  finite-difference  method 
The  original  problem  is 

F(x)  =  I(x)  +  VQ(x)  =  0  (4.9) 

where  x  is  the  circuit  waveform  vector,  I  is  the  contribution  from  the  resistive  part, 
Q  is  the  contribution  from  the  reactive  part  andV  represents  the  discretized  time 
derivative.  The  homotopy  technique  sweeps  from  an  initial  value  problem 
(A  =  0)  to  a  periodic  boundary  value  problem  (A  =  1) : 


F(x,Ax,  co)  =  I(x)  +  V<2(.v)  +  V<2(.v)  =  0 
x  =  \x  +  (1  -  A,  )x0 


(4.10) 


where  V  and  V  are  block  lower  triangular  and  strictly  block  upper  triangular 
matrices,  respectively,  and  V  =  V  +  V  .  This  homotopy  mapping  has  been  applied 
to  solve  different  oscillators  where  Newton’s  method  fails. 


41 


4.4  Homotopy  Method  Solver:  HOMPACK 

HOMPACK  is  a  large  collection  of  Fortran  subroutines  for  solving  various 
nonlinear  systems  of  equations  by  homotopy  methods  [11].  The  algorithms  in 
HOMPACK  are  based  on  “artificial-parameter  generic  homotopies”.  While 
Newton’s  method  only  has  a  local  convergence,  HOMPACK  is  globally 
convergent  with  probability  one  under  mild  assumptions.  These  assumptions  are 

-  A  unique  solution  at  >.=0. 

-  A  bifurcation  free  zero  curve. 

-  A  bounded  solution  without  rapid  oscillation  along  the  zero  curve. 

At  the  end  of  zero  curve  (when  h=  1 )  singularities  may  be  encountered.  This 
happens  precisely  when  the  original  problem  is  singular  at  the  solution.  For  some 
mild  singularities,  the  homotopy  method  can  remain  nonsingular  at  the  solution, 
but  in  the  general  case  this  is  not  so. 

There  are  three  different  algorithms  in  HOMPACK: 

-  Fixpdf:  ordinary  differential  equation-based  method. 

-  Fixpnf:  augmented  Jacobian  method. 

-  Fixpqf:  normal  flow  algorithm. 

Among  them,  the  fixpdf  algorithm  is  the  most  robust  but  it  generally  takes  more  k 
steps  and  so  it  is  computationally  more  expensive.  It  is  difficult  to  predict  which 
algorithm  will  be  the  best  for  a  given  problem. 

HOMPACK  divides  the  problem  into  three  subdivisions: 

■  Fixed-point  problem 

The  equation  x  =  Fix)  is  solved  by  following  the  zero  curve  of  the 


homotopy  map 
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H(x,  A)  =  (1  -  A)(x  -a)  +  A(x  -  F(x))  (4. 1 1) 

starting  from  A  =  0  and  x  =  a  .  The  curve  is  parameterized  by  arc  length  5. 

■  Zero-finding  problem 

The  equation  F(x)  =  0  is  solved  by  following  the  zero  curve  of  the 
homotopy  map 

H(x,A)  =  (l-A)(x-a)  +  AF(x )  (4.12) 

starting  from  A  =  0  and  x  =  a  .  The  curve  is  also  parameterized  by  arc  length  s. 


■  Curve-tracking  problem 

The  curve-tracking  problem  deals  with  an  arbitrary  embedding  of  A.  The 
homotopy  map  H(a,A,x)  is  assumed  to  satisfy 

dH(a,A,x )  dH(a,A,x ) 

rank[ - — - , - - - ]  =  N  (4.13) 

d/t  dx 

for  all  points  ( A ,  x)  such  that  H(a ,  A,x)  =  0 .  It  is  further  assumed  that 

mnk[dH(a,Xx)]  =  N  (414) 

dx 

with  a  fixed.  The  zero  curve  of  H(a,A,x )  starting  from  A  =  0  and  x  =  x0  is 
tracked  until  A  =1  by  solving  the  ordinary  differential  equation 
d  H  ( a,A(s),x(s )) 


ds 


=  0  where  s  is  the  arc  length  along  the  zero  curve.  Also  the 


homotopy  map  FI  (a ,  A,  x)  =  0  is  assumed  to  be  constructed  such  that  ^  =  0  . 

ds 
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For  the  fixed-point  and  zero-finding  problems,  the  user  must  supply  a 
subroutine  f(x,v )  which  evaluates  Fix)  at  x  and  returns  the  value  to  vector  v, 
and  a  subroutine  fjac(x,  v,k)  which  returns  in  v  the  Ath  column  of  the  Jacobian 
matrix  of  F(x )  evaluated  at  x. 


For  the  curve  tracking  problem,  the  user  must  supply  a  subroutine 
rho(a,A,x,v)  which  evaluates  H(a,; l,x)  at  (a,/ l,x)  and  returns  the  value  to 
vector  v,  a  subroutine  rhoa(a,A,x,  par,ipar)  which  given  (A,  x)  returns  a 
parameter  vector  in  a  such  that  Hia.A.x)  =  0.  and  a  subroutine 
rhojacia,  A,  x,v,k,  par,  ipar)  which  returns  in  v  the  Ath  column  of  the  N  x  (N  + 1) 
dH(a,A,x )  dH(a,A,x ) 


Jacobian  matrix  [- 


dA 


dx 


-]  evaluated  at  ( a,A,x ). 


4.5  Homotopy  Formulations  for  Frequency-Domain  Simulation  of 

Oscillators 

As  has  been  shown  in  Section  4.1,  obtaining  an  initial  guess  for  the  probe 
voltage  for  some  high-Q  oscillators  is  difficult.  For  this  reason,  Newton’s  method 
can’t  be  used  to  solve  the  top-level  probe  equations.  The  application  of  a 
homotopy  method  is  described  in  this  section. 

The  original  equations  are  the  top-level  probe  equations  as  shown  in  Eq. 
(3.10).  In  order  to  apply  a  homotopy  map,  Eq.  (3.10)  first  needs  to  be  modified  by 
dividing  with  the  probe  voltage  as  shown  in  Eq.  (4.15).  The  purpose  of  doing  so  is 
to  isolate  the  real  solution  from  the  trivial  DC  solution  of  Eq.  (3.10)  which  the 
homotopy  map  is  more  easily  attracted  to. 
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Kro^pro.,'’^ 

v 

probe 

isProbe(yprobe,co) 


V 


=  0 


=  0 


probe 


(4.15) 


A  homotopy  map  is  then  applied  to  the  above  modified  top-level  probe 
equations. 

V  .  (V  ,  ,a) 

x  pr»^  robe  +{l_A)gi  (vu  _  )  =  0 

probe  (4.16) 

7s  .  (V  .  , (O) 

*  P  yP  +  (1  -  A)g2(a>-  a>inii )  =  0 

probe 


where  Vm,t  and  (!),na  are  the  initial  values  for  the  probe  voltage  and  frequency,  and 
gi  and  g2  are  scaling  factors.  The  typical  value  for  g/  is  10  4  and  y?  is  chosen  such 
that  both  terms  in  Eq.  (4.16)  are  of  a  similar  magnitude  depending  on  the 
frequency  of  oscillation.  From  Eq.  (4.16)  we  can  see  that  by  sweeping  from  A  =  0 
to  A  =  1 ,  the  homotopy  map  goes  gradually  from  an  easy  to  solve  forced  circuit 
problem  to  a  hard  to  solve  problem,  i.e.,  the  original  oscillator  problem. 

A  Pierce  crystal  oscillator  from  [8,  9]  is  shown  in  Fig.  4.2  and  used  to 
evaluate  the  above  homotopy  algorithm.  The  crystal  Q  is  2.5xl04.  The  circuit  has 
two  poles  (1.6774645xl04  ±  j6.289348xl07),  and  the  initial  guess  for  the 
oscillation  frequency  is  set  to  10.008745MHz  correspondingly.  The  voltage  probe 
is  placed  at  the  high- impedance  collector  node  of  the  BJT.  Figs.  4.3  and  4.4  show 
the  solution  traces  with  the  initial  probe  voltage  starting  from  a  probe  voltage 
equal  to  the  DC  operating  point  value  and  lmV,  respectively.  The  traces  are  from 
the  FIXPDF  algorithm.  Similar  traces  are  obtained  with  the  FIXPQF  and  FIXPNF 
algorithms. 
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Figure  4.2:  Circuit  schematic  of  the  Pierce  oscillator.  R1  =  100KQ.  R2=2. 
Cl=100pF,  C2=100pF,  Cp=25pF,  Cc=99.5fF,  Rc=6.4Q,  Lc=2.5 

Vcc=12V,  p  =100  for  the  BJT. 


Figure  4.3:  Solution  traces  for  the  Pierce  oscillator  with  the  initial  probe 
amplitude  equal  to  the  DC  solution,  (a)  Probe  voltage  as  a  function  of  X.  (b) 
Oscillation  frequency  as  a  function  of  X. 


<N  in 
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(a) 


Figure  4.4:  Solution  traces  for  the  Pierce  oscillator  with  the  initial  probe 
amplitude  equal  to  lmV.  (a)  Probe  voltage  as  a  function  of  X.  (b)  Oscillation 
frequency  as  a  function  of  X. 


From  Figs.  4.3  and  4.4  we  can  see  the  progress  of  the  probe  voltage  and 
oscillation  frequency  as  a  function  of  the  embedding  parameter  X.  For  the  curve 
starting  from  the  DC  operating  point  value,  the  probe  voltage  tends  to  go  to  the 
DC  solution  (a  probe  amplitude  of  OV),  but  the  modified  equations  handle  this 
situation  and  the  probe  voltage  converges  to  a  negative  value  that  is  a  valid 
solution.  For  the  curve  starting  from  a  small  value,  the  tracking  curve  goes  to  the 
periodic  steady-state  solution  directly  since  the  original  probe  equations  have 
been  effectively  modified  to  isolate  the  trivial  DC  solution. 

The  problem  with  the  embedding  in  Eq.  (4.  16)  is  that  it  converges  slowly 
although  it  is  robust.  In  this  particular  case,  more  than  a  thousand  steps  (several 
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thousand  Jacobian  evaluations)  were  needed  to  obtain  the  solution.  A  closer 
examination  of  the  solution  traces  reveals  the  possible  causes.  The  probe  voltage 
and  frequency  are  two  different  parameters  and  embedding  them  with  the  same 
parameter  X  is  not  the  most  efficient  approach.  Furthermore,  from  Fig.  4.3  we  see 
that  from  X  =  0  to  about  X  =  0.7,  the  change  in  frequency  is  very  small,  but  the 
probe  voltage  progresses  quickly  even  with  turning  points.  After  X  =  0.7,  the 
solution  trace  for  frequency  shows  a  change  and  the  tracking  takes  several  steps. 
Fig.  4.4  shows  the  same  trend.  This  suggests  that  the  tracking  for  the  voltage  and 
the  oscillation  frequency  should  be  decoupled  by  using  two  different  embedding 
parameters.  For  this  reason,  ki  is  used  for  embedding  the  probe  voltage  and  X2  is 
used  for  embedding  the  oscillation  frequency.  Then  a  revised  algorithm  is 
obtained  as  shown  in  Fig.  4.5. 


While  (X#l)  do  { 

1.  For  each  X\  (fixed  probe  voltage),  track  X2  from  0  to  1. 

7s  .  (V  ,  ,CO) 

K  P  P  +(1  -*2)g2{co- (Oinit)  =  0 

probe 

2.  When  X2  equals  1  (fixed  frequency),  step  X\  . 

V  .  (V  ,  ,0)) 

\  Pr0t ^  ^  +  a  -  4  )gl  (Vpmbe  ~  )  =  0 

probe 

} 

Solution  is  obtained  for  Xi=l. 


Figure  4.5:  Algorithm  1:  A  two-parameter  embedding  homotopy  algorithm. 
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Step  w 

(b) 

Figure  4.6:  The  solution  tracking  process  in  a  high-Q  resonator,  (a)  The  voltage 
in  the  outer  loop  and  frequency  in  the  inner  loop,  (b)  The  frequency  in  the  outer 
loop  and  voltage  in  the  inner  loop. 

In  this  algorithm  X\  (probe  voltage)  and  Ao  (frequency)  are  solved  iteratively. 
A]  (probe  voltage)  is  tracked  in  the  outer  loop  while  Ai_  (frequency)  is  tracked  in 
the  inner  loop.  The  reason  for  this  ordering  can  be  see  from  the  tracking  process 
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in  a  high-Q  resonator  as  shown  in  Fig.  4.6.  From  Fig.  4.6(a)  we  can  see  if  the 
voltage  is  in  the  outer  loop,  at  a  fixed  voltage  Vi  we  will  get  the  frequency 
solution  (O j,  then  the  voltage  is  stepped  to  a  new  value  V2  using  the  embedding 
function  as  in  Fig.  4.5.  For  this  new  voltage  V2,  we  need  to  solve  for  the 
corresponding  co2  .  As  seen  from  Fig.  4.6(a),  the  frequency  corresponding  to 
voltage  V 1  serves  as  a  good  initial  guess  since  it  is  quite  close  to  co2  .  However,  if 
the  process  is  reversed  as  in  Fig.  4.6(b)  where  the  frequency  is  in  the  outer  loop 
and  the  voltage  is  in  the  inner  loop,  a  very  small  step  in  a  (from  CO]  to  co2 )  will 
introduce  a  large  change  in  the  voltage  (from  Vi  to  V2).  In  this  case,  Vi  is  no 
longer  a  good  initial  guess  for  the  new  solution  V2\ 

The  globally  convergent  algorithms  in  [11]  are  based  on  a  single  homotopy 
parameter  X.  The  above  algorithm  is  used  with  HOMPACK  by  switching  the  two 
embedding  parameters.  This  algorithm  provides  a  significant  speed  up  even 
though  A]  and  X2  are  solved  iteratively.  This  algorithm  can  be  further  modified  by 
noting  that  for  a  small  value  of  the  initial  probe  voltage  (for  starting  the  homotopy 
curve  tracking),  the  oscillation  frequency  estimated  by  a  linear  pole-zero  analysis 
serves  as  a  very  good  initial  guess.  Therefore,  a  homotopy  mapping  to  track  the 
frequency  is  not  necessary  and  one  can  simply  use  a  Newton  method  in  Step  1  of 
Algorithm  1.  By  doing  so,  the  modified  algorithm  can  be  made  more  efficient 
since  the  homtopy  method  is  generally  more  expensive  than  the  Newton  method. 
In  addition,  only  one  embedding  parameter  is  needed  which  makes  the 
implementation  much  easier.  The  corresponding  algorithm.  Algorithm  2,  is  given 
in  Fig.  4.7. 


1  For  high  Q  resonators,  the  voltage  changes  significantly  for  a  small  shift  in  frequency. 
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The  new  algorithms  have  been  tested  on  different  high-Q  and  moderate  Q 
oscillators  to  demonstrate  their  convergence  performance.  The  results  are 
summarized  in  Chapter  7. 


The  probe  voltage  starts  from  a  small  value. 

While  (A^l)  do  { 

3.  For  each  X  (fixed  probe  voltage),  solve 

=  0 

4.  Step  X 

r  av  .,cd) 

*  '  y  '  +  (1  ~  ^  )ffl  (Y probe  ~  Vini,  )  =  0 

probe 

} 

Solution  is  obtained  for  A^l. 


Figure  4.7:  Algorithm  2:  Modified  embedding  algorithm  for  high-Q  oscillator 
simulation. 


4.6  Summary 

In  this  chapter  globally  convergent  homotopy  methods  have  been  applied  to 
address  the  frequency-domain  simulation  of  high-Q  oscillators.  The  basic  theory 
of  homotopy  methods  and  its  application  to  DC  operating  points  and  time-domain 
steady-state  oscillator  simulation  are  reviewed.  A  brief  description  of  the 
homotopy  method  solver,  HOMPACK,  is  also  given.  Different  embedding 
techniques  have  been  proposed  for  frequency-domain  high-Q  oscillator  simulation 
and  an  efficient  and  robust  algorithm  is  identified. 
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5.  FREQUENCY-DOMAIN  SIMULATION  OF  HIGHLY  NONLINEAR 

RING  OSCILLATORS 


5.1  Introduction 

All  previous  work  using  the  harmonic  balance  method  for  oscillator 
simulation  has  focused  on  LC  oscillators  which  have  near  sinusoidal  waveforms 
[7-10,  19].  Some  of  these  techniques  can  be  applied  to  simulate  ring  oscillators, 
but  due  to  the  highly  nonlinear  nature  of  these  circuits,  the  simulation  fails  to 
converge.  The  mechanism  of  divergence  for  ring  oscillators  is  different  from  that 
of  high-Q  oscillators.  Ring  oscillators  have  a  low  Q  and  the  solution  of  highly 
nonlinear  circuit  equations  is  a  problem  rather  than  the  initial  guess  for  the  probe 
frequency  and  voltage. 

Commercially  available  harmonic  balance  simulators  can  only  simulate 
highly  nonlinear  ring  oscillators  using  results  from  a  time-domain  transient 
analysis  as  an  initial  guess.  This  approach,  however,  constrains  the  simulation  of 
ring  oscillators  in  a  RF  system  where  transient  analysis  is  no  longer  efficient 
when  mixer  circuits  are  present.  In  this  chapter,  two  new  algorithms  for  robust 
simulation  of  ring  oscillators  are  proposed. 

This  chapter  is  organized  as  follows.  The  challenges  in  ring  oscillator 
simulation  in  the  frequency  domain  are  introduced  in  Section  5.2.  A  phase 
relationship  for  ring  oscillators  is  identified  in  Section  5.3.  By  taking  advantage  of 
the  ring  oscillator  phase  relationship,  two  new  methods  for  frequency-domain  ring 
oscillator  simulation  are  proposed.  A  single-delay  cell  method  is  described  in 
Section  5.4  and  a  multiple-probe  method  is  discussed  in  Section  5.5.  This  chapter 
is  summarized  in  Section  5.6. 
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5.2  Convergence  of  Harmonic  Balance  Methods  for  Ring  Oscillators 

A  comparison  is  made  between  the  performance  of  the  traditional  harmonic 
balance  method  (HB)  in  Eq.  (3.6),  the  continuation  method  of  [8]  (Conti),  the 
modified  continuation  method  of  [9]  (Cont2)  and  the  two-level  method  of  [7]  for 
different  stage  ring  oscillators  with  a  basic  inverter  cell  (Fig.  5.1  (a)).  The  BSIM3 
MOSFET  model  is  used  for  all  MOS  devices.  The  traditional  harmonic  balance 
method  and  continuation  method  1  (Conti)  failed  after  300  iterations  for  all  cases. 
The  modified  continuation  method  (Cont2)  worked  only  for  the  first  case. 
However,  it  is  computationally  much  more  expensive  and  less  robust  than  the 
two-level  method  for  ring  oscillator  simulation. 


Table  5.1:  Convergence  performance  comparison  of  different  harmonic  balance 
methods  for  ring  oscillator  simulation,  x  denotes  no  convergence  in  300 
iterations. 


#  stage 

#Freq  for 
HB 

Analysis 

Freq  (GHz) 

#  Iterations 

HB 

Conti 

Cont2 

2-level 

3 

16 

3.25 

X 

X 

187 

23 

5 

32 

1.95 

X 

X 

X 

38 

7 

64 

1.4 

X 

X 

X 

43 

Although  the  two-level  method  is  more  robust  than  the  continuation  method 
for  ring  oscillator  simulation,  its  convergence  is  still  very  sensitive  to  the  initial 
guess.  For  some  circuits  the  initial  frequency  guess  must  be  within  1%  of  the 
oscillation  frequency  to  reach  convergence.  Other  circuits  do  not  converge  even 
when  the  exact  oscillation  frequency  and  probe  voltage  are  given  as  the  initial 
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guess.  This  is  illustrated  with  results  summarized  in  Table  5.2.  In  this  table, 
\Af\/ fosc  ’\AVProbe\/Vosc  are  the  errors  of  the  initial  guess  relative  to  the  solution  for 

the  frequency  and  probe  voltage,  respectively. 


Table  5.2:  Convergence  properties  of  different  ring  oscillator  circuits,  x  denotes 
no  convergence  under  exact  initial  guess  for  oscillation  frequency  and  probe 
voltage. 


Oscillator  circuit 

Circuit  1 

Circuit  2 

Circuit  3 

Circuit  4 

Delay  cell  (Fig.  5.1) 

(a) 

(b) 

(c) 

(d) 

#  Stages 

9 

21 

3 

5 

|A/| ffosc 

X 

1% 

X 

1% 

\wprobe\/vosc 

X 

5% 

X 

5% 

A  close  examination  of  the  results,  for  ring  oscillators  for  which  the 
simulation  always  fail,  reveals  convergence  problems  in  the  solution  of  the 
bottom-level  voltage-probe  forced  circuit.  A  circuit  waveform  (the  cosine  part  of 
the  first  harmonic  of  a  node  voltage)  during  the  Newton  iterative  process  from 
Circuit  3  is  shown  in  Fig.  5.2.  From  Fig.  5.2  (a)  we  can  see  that  the  circuit 
waveform  diverges  even  with  a  very  small  damping  factor  for  the  Newton  method. 
To  overcome  this  problem,  voltage  and  current  thresholds  can  be  set,  but  from  Fig. 
5.2  (b)  we  see  that  the  circuit  waveform  oscillates  between  these  thresholds  and 
fails  to  converge. 
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VDD  VDD 


Figure  5.1:  Ring  oscillator  delay  cells,  (a)  Basic  delay  cell,  (b)  Current  starved 
delay  cell .  (c)  Vanilla  delay  cell,  (d)  Cross-coupled  load  delay  cell  [39]. 
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(a) 


(b) 


Figure  5.2:  Circuit  waveform  during  the  Newton  iterative  process,  (a)  Use  of  a 
small  damping  factor,  (b)  Adaptive  damping  with  voltage/current  thresholds. 
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The  use  of  a  homotopy  method  at  the  bottom-level  has  been  investigated. 
The  divergence  in  the  ring  oscillator  simulation  comes  from  the  nonlinear  active 
devices.  So  an  embedding  function  that  starts  from  only  linear  capacitances  or 
linear  capacitances  and  resistances  has  been  tried.  However,  this  approach  did  not 
give  a  good  performance. 

When  a  voltage  source  with  the  exact  oscillation  waveform  at  a  node  is  used 
to  force  the  ring  oscillator  circuit  at  that  node,  the  circuit  converges  readily.  This 
suggests  that  besides  the  high  nonlinearity  of  the  circuit,  the  manner  in  which  the 
ring  oscillator  circuit  is  converted  to  a  forced  circuit  also  plays  a  critical  role  in 
the  convergence  performance. 

Based  on  the  above  observations  two  new  approaches  are  proposed  for 
robust  ring  oscillator  simulation.  The  first  approach  relies  on  a  decrease  in  the 
number  of  nonlinear  devices.  The  second  approach  uses  a  more  effective  forcing 
function.  Both  approaches  exploit  the  periodic  structure  of  ring  oscillators. 


5.3  Ring  Oscillator  Phase  Relationships 

Ring  oscillators  have  a  periodic  structure,  whereby  the  waveforms  at  the 
different  nodes  differ  in  phase  by  multiples  of  Infm ,  where  m  is  the  number  of 
the  stages  in  the  oscillator  [40].  The  explicit  phase  shift  relationship  between  the 
input  and  output  nodes  of  each  delay  cell  can  be  seen  from  the  time  domain 
waveforms  of  a  five-stage  ring  oscillator  as  shown  in  Fig.  5.3.  In  this  particular 

case,  the  input  and  output  differ  in  phase  by  ^-(216°). 
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Figure  5.3:  Five- stage  ring  oscillator  waveform. 


A  general  relationship  for  the  phase  shift  of  each  delay  cell  is  now  derived. 
Recall  that  the  period  of  oscillation  for  a  ring  oscillator  is  given  by  [35,  40] 

T  =  2mTD  (5.1) 

where  rD  is  the  time  delay  of  the  delay  cell.  From  Eq.  (5.1),  the  time  delay  rD 
corresponds  to  a  phase  shift  of  k /m  radians.  Since  the  inverting  delay  cell 
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contributes  an  additional  phase  shift  of  k  ,  the  phase  delay  for  each  inverting  cell 
can  be  expressed  as 

e  =  (— )^  (5.2) 

m 

The  phase  relationship  for  a  three-stage  ring  oscillator  is  shown  in  Fig.  5.4  which 
pictorially  shows  the  above  derivation.  From  Eq.  (5.2)  the  total  phase  shift  of  a  m 
stage  ring  oscillator  will  be  {m+X)7t  when  m  is  an  odd  number.  If  m  is  an  even 
number,  and  noting  that  the  last  stage  of  the  ring  oscillator  is  non-inverting,  the 
total  phase  shift  is  mn  instead.  In  both  cases,  the  total  phase  shift  is  a  multiple  of 
2 n  which  guarantees  a  positive  feedback. 


Figure  5.4:  Phase  relationship  for  a  three-stage  ring  oscillator. 
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Based  on  the  above  observation  we  have  the  input-output  voltage  (current) 
phase  relationship  of  a  single-ended  ring  oscillator  delay  cell  as 


-jM 


(5.3) 


where  X  is  the  frequency-domain  circuit  waveform  (complex  vector),  k  is  the 
harmonic  order,  and  6  is  the  phase  relationship  given  in  Eq.  (5.2).  For  differential 
ring  oscillators,  this  relationship  holds  for  each  half  circuit,  while  the  phase 
difference  between  the  two  half  circuits  is  simply  jt. 


5.4  Single-delay  Cell  Method 

The  single-delay  cell  method  for  ring  oscillator  simulation  using  the 
harmonic  balance  method  recognizes  the  structure  of  the  oscillator  assuming  all 
cells  are  identical.  The  phase  relationships  between  the  individual  delay  cells 
derived  in  the  previous  section  are  used  to  simulate  a  single-delay  cell  and 
determine  the  periodic  steady  state  of  the  oscillator. 

5.4.1  Single-delay  Cell  Equivalent  Circuit 

The  basic  idea  behind  the  single-delay  cell  method  is  that  one  delay  cell  as 
well  as  the  information  about  the  number  of  stages  gives  enough  information 
about  the  operation  of  the  ring  oscillator  circuit.  Therefore,  instead  of  simulating 
the  complete  ring  oscillator  circuit,  we  can  simplify  the  simulation  to  that  of  a 
single-delay  cell  equivalent  circuit  as  shown  in  Fig.  5.5.  The  delay  element  gives 
the  phase  relationship  between  the  input-output  waveforms  according  to  Eq.  (5.3), 
which  is  a  function  of  the  number  of  stages  as  shown  in  Eq.  (5.2). 
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Figure  5.5:  Single-delay  cell  equivalent  circuit  for  ring  oscillator  simulation. 


For  differential  ring  oscillators,  this  relationship  holds  for  each  half  circuit, 
while  the  phase  difference  between  the  two  half  circuits  is  simply  n.  The 
harmonic  balance  method,  being  a  frequency-domain  method,  can  handle  the 
phase  delay  elements  easily  and  efficiently.  With  a  reduction  in  the  number  of 
nodes  in  the  circuit,  the  simulation  is  considerably  sped  up.  Furthermore,  the 
convergence  is  also  improved  since  fewer  nonlinear  devices  are  included  in  the 
simulation. 


5.4.2  Matrix  Formulation  for  Harmonic  Balance  Method 

The  input  and  output  voltage/current  complex  vectors  can  be  expressed  in 
rectangular  coordinates  as 


X\  =  xkt 

out  out,c 


I  •  K 

+  JXo, 


Xk=xk  +  jxk 

in  in.c  J  u 


(5.4) 


Substituting  Eq.  (5.4)  into  Eq.  (5.3)  we  see  that 
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•Ce  =  cos  (kd)xkinc  +  sin  {kd)xkins  ^ 

xk ,  = -sin(/c#);C  +cos(/c^)xA: 

where  the  subscripts  c  and  .s'  denote  the  cosine  and  sine  parts  of  the  circuit 
waveform  Fourier  coefficients,  respectively.  The  voltage  (current)  phase  delay 
can  be  stamped  similar  to  the  branch  constitutive  equations  directly  in  the 
harmonic  balance  Jacobian  matrix  as  shown  in  Eq.  (5.6) 


Nin 

Nqu  t  I in 

Nin  Nout 

Iin 

Nin  Noul 

h 

N, 

X 

x  -lj 

DC 

N  , 

out 

X 

x  i  | 

1 

branch 

-1 

l  ! 

_i 

Stamp  = 

(5.6) 

N. 

m 

X  X 

-l  ! 

COSk 

Nmt 

X  X 

cos  kO  | 

i 

sin  kd 

branch 

-  cos  kO  1 

i 

i 

-  sin  kd 

Nin 

X  X 

-1  j 

SINk 

-  sin  kO 

X  X 

cos  kd 

branch 

sin  kd 

-  cos  kd  1 

1 

1 

1 

The  branch  constitutive  relations  (BCR)  specify  the  constraints  of  the  “delay(k60” 
branch  in  Fig.  5.5.  With  this  Jacobian  matrix  we  can  solve  the  ring  oscillator 
single-delay  cell  equivalent  circuit  using  Eq.  (3.6).  However,  a  straightforward 
implementation  suffers  from  convergence  problems  since  the  frequency  is  an 
unknown.  To  overcome  the  convergence  problem  we  use  the  voltage  probe 
technique  from  [7]. 
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5.4.3  Modified  Voltage  Probe  Method 

The  voltage  probe  is  used  to  convert  the  single-delay  cell  into  a  forced 
equivalent  circuit.  At  f=fOSc,  the  input  of  the  delay  cell  is  driven  with  a  voltage 
probe,  and  the  input-output  current  still  satisfy  the  relationship  as  in  Eq.  (5.3). 
This  is  expressed  as  a  CCCS  in  the  circuit  schematic  of  Fig.  5.6.  For  frequencies 
fffosc,  the  circuit  is  the  same  as  that  in  Fig.  5.5.  The  solution  is  obtained  when  the 
input-output  voltages  at  the  fundamental  frequency  satisfy  Eq.  (5.3). 


Figure  5.6:  Single-delay  cell  equivalent  circuit  with  a  voltage  probe  at  f  =  fosc . 


5.4.4  Solution  by  Newton’s  Method 

The  circuit  is  now  solved  with  a  two-level  method.  The  first  level  of  the 
two-level  method  solves  the  harmonic  balance  equations  of  the  voltage  probe 
forced  circuit. 


'  K 

-<d) 

'x°'+1)l 

RhsF  U) " 

«X)T 

0 

0 

AVC 

= 

v 

probe 

<(Dr 

0 

0 

1 

> 

1 _ 

0 

(5.7) 
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Ay  =-cos^y1  -  sin  dVl  +  V1 

c  in,c  in,s  ou 

AV  =  sin  0Vl  -costfE1  +V\ 

s  in,s  in,c  out,s 


(5.8) 


where  j u  is  the  modified  harmonic  balance  Jacobian  with  the  branch  constitutive 
equations  stamped  according  to  Eq.  (5.6).  e'in(\)  (A  (1))  is  the  unit  vector  that 
selects  the  cosine  (sine)  part  of  the  first  harmonic  of  the  node  to  which  the 
probe  is  attached.  ec.  (1)  (<?)  (1))  is  the  unit  vector  that  selects  the  branch 
constitutive  equations  in  Eq.  (5.6)  to  open  the  voltage  delay  circuit  at  the 
fundamental  frequency.  V'n  c  ( V'n  s )  is  the  cosine  (sine)  part  of  the  first  harmonic  of 

the  input  voltage.  Vo1utc  <  V'mt  5 )  is  the  cosine  (sine)  part  of  the  first  harmonic  of  the 
output  voltage. 

The  second  level  solves  the  probe  equations: 


AK(Vpro^)  =  0 


(5.9) 


The  Jacobian  for  this  set  of  equations  is: 


J 


probe 


d\vc(yprobe,co)  d\vc{vprobe,co) 

dVPro„e  C)0J 

dAVs(Vpm„e,a)  dAVs(Vprobe,C0) 

dVprobe  dco 


(5.10) 


This  2x2  Jacobian  can  be  obtained  by  a  simple  computation  using  information 
from  the  first-level  solution. 

The  frequency-domain  nonlinear  equations  for  the  circuit-probe 
combination,  G,  can  be  expressed  as: 

G(Xa,cd)  =  0 


(5.11) 
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where  xA  is  the  vector  of  the  circuit  unknowns,  and  AVC  and  AVS.  The  sensitivity 
with  respect  to  Vprobe  and  co  can  then  be  computed  as: 


3 G  dXA  _  3 G 

a g  dxA  __a g 

dX  .  dco  d&) 

A 


(5.12) 


where  - is  the  augmented  Jacobian  in  Eq.  (5.7).  All  entries  of  the  2x2  Jacobian 

dXA 

matrix  in  Eq.  (5.10)  can  be  extracted  from  the  solution  of  Eq.  (5.12). 


5.4.5  Differential  Ring  Oscillator  Simulation 

For  differential  ring  oscillators,  the  voltage  probe  is  added  to  the  circuit  at 
the  oscillation  frequency  differentially.  By  doing  this,  the  phase  constraints  are 
automatically  forced  for  the  two  half  circuits.  Each  half  circuit  has  the  phase 
relationship  as  in  the  single-ended  oscillator.  The  equivalent  schematic  of  the 
circuit  in  the  frequency  domain  is  shown  in  Figs.  5.7  and  5.8. 

To  effectively  open  the  voltage  delay  block  for  both  half  circuits,  we 
observe  that  AV  and  AV  for  the  two  half  circuits  have  a  phase  difference  of  jt. 
Using  this  information,  the  Jacobian  can  be  constructed  as: 

"  JM  -/;(!)  -//(l)Txa+1)l  \RhsFur 
W)r  0  o  AV,  =  (5.13) 

Is  (l)r  0  0  AV  0 

where  If  1)  (/:’(!))  is  of  the  form  [...  1  ...  -1  ...f  ,  which  selects  the  branch 
constitutive  equations  in  Eq.  (5.6)  to  open  the  voltage  delay  circuit  at  the 


65 


fundamental  frequency  for  each  half  circuit.  /‘n(lj  (/'(l))  is  also  of  the  form 
[...  1  ...  -1  ...f  where  the  cosine  (sine)  part  of  the  first  harmonic  of  the 
differential  input  nodes  to  which  the  probe  is  attached  are  selected. 


Figure  5.7:  Single-delay  cell  equivalent  circuit  for  a  differential  ring  oscillator  at 

ftfosc. 


Figure  5.8:  Single-delay  cell  equivalent  circuit  with  a  voltage  probe  for  the 
differential  ring  oscillator  at  f=fosc. 


66 


5.4.6  Limitation 

The  proposed  single-delay  cell  method  assumes  that  all  the  oscillator  cells 
are  identical  and  cannot  be  used  when  large  mismatches  between  cells  need  to  be 
accounted  for.  More  work  needs  to  be  done  on  extensions  of  this  work  to  non¬ 
identical  delay  cells. 


5.5  Multiple-probe  Method 

The  multiple-probe  method  for  ring  oscillator  simulation  using  the  harmonic 
balance  method  recognizes  that  the  manner  in  which  the  ring  oscillator  circuit  is 
converted  to  a  forced  circuit  plays  a  critical  role  in  convergence.  A  more  effective 
way  of  converting  the  ring  oscillator  circuit  to  a  forced  circuit  is  accomplished 
with  the  help  of  multiple  probes.  Although  multiple  probes  are  applied  in  the 
circuit,  the  initial  guess  still  remains  the  same  as  that  of  the  single-probe  method 
by  taking  advantage  of  the  ring  oscillator  phase  relationships. 


5.5.1  Application  of  Multiple  Probes 

To  show  how  the  divergence  in  the  bottom- level  probe  forced  circuit  occurs, 
the  Newton  iterative  process  for  the  fundamental  and  2nd  harmonic  at  different 
nodes  for  a  single-ended  9-stage  ring  oscillator  (Circuit  1  in  Table  5.2)  is  shown 
in  Figure  5.9.  The  Newton  iterative  process  starts  from  some  initial  value.  At  the 
first  iteration  the  fundamental  voltages  at  some  nodes  go  to  a  large  number  while 
all  the  2nd  harmonics  are  still  reasonable.  However,  at  the  second  iteration  the  2nd 
harmonics  diverge  quickly  due  to  the  divergence  of  the  fundamental  circuit 
waveform  at  the  previous  iteration.  The  same  trend  can  be  seen  with  other  higher 
order  harmonics. 
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Figure  5.9:  Newton  iterative  process  for  the  fundamental  and  2nd  harmonic  at 
different  nodes  and  different  iterations,  (a)  initial  guess,  (b)  iteration  #1,  (c) 
iteration  #2. 


Based  on  above  observation,  a  new  way  of  forcing  the  ring  oscillator  circuit 
has  been  developed  using  multiple  probes.  A  voltage  probe  is  applied  at  the  input 
of  each  delay  cell  as  shown  in  Fig.  5.10.  By  doing  so,  the  fundamental  voltage 
waveform  at  all  the  important  nodes  in  ring  oscillator  circuits  are  forced,  which 
helps  in  constraining  the  fundamental  waveform  and  thus  avoiding  convergence 
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problems  as  demonstrated  in  Fig.  5.9.  For  differential  ring  oscillators,  voltage 
probes  are  added  to  the  circuit  at  the  oscillation  frequency  differentially. 


Figure  5.10:  Ring  oscillator  simulation  using  multiple  voltage  probes. 


For  most  ring  oscillator  circuits,  there  is  no  need  to  apply  as  many  voltage 
probes  as  the  number  of  stages.  Simulations  that  fail  with  a  single  probe  can 
converge  readily  with  the  application  of  two  or  more  voltage  probes.  However, 
some  ring  oscillator  circuits  do  require  m  (the  number  of  the  stages  in  the  ring 
oscillator)  voltage  probes  for  reliable  convergence.  For  a  certain  ring  oscillator 
structure  it  is  difficult  to  determine  the  minimum  number  of  voltage  probes  or 
where  these  probes  should  be  placed.  For  this  reason  voltage  probes  are  generally 
applied  at  the  input  of  each  stage.  Although  a  larger  number  of  probes  are  used, 
the  increase  in  the  computational  effort  at  each  bottom  level  iteration  is  small  as 
shown  later  in  Section  5.5.4. 
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5.5.2  Initial  Guesses  for  Multiple  Probes 

Since  each  voltage  probe  requires  its  own  initial  guess,  multiple  probes 
require  multiple  guesses,  which  can  be  an  impractical  task  for  most  circuits. 
However,  ring  oscillators  have  a  periodic  structure  and  when  all  delay  cells  are 
identical  the  waveforms  at  the  different  nodes  differ  only  in  phase  by  an  angle 
given  in  Eq.  (5.2).  If  one  voltage  probe  is  chosen  as  a  reference,  the  initial  guess 
for  all  other  voltage  probes  can  be  obtained  by  a  phase  shift  using  Eq.  (5.3).  Even 
when  there  are  mismatches  between  cells,  Eq.  (5.3)  still  serves  as  a  good 
approximation.  For  differential  ring  oscillators,  the  phase  relationship  in  Eq.  (5.3) 
holds  for  differential  waveforms. 


5.5.3  Solution  by  Newton’s  Method 

The  multiple-probe  method  is  also  implemented  by  a  two-level  Newton 
procedure.  The  bottom  level  of  the  two-level  method  solves  the  harmonic  balance 
equations  of  the  multiple-probe  forced  circuit.  The  sine  part  of  one  probe  voltage 
is  set  to  zero  to  fix  the  phase  and  this  probe  voltage  is  also  used  as  a  reference  for 
the  initial  guess  of  the  other  probes.  The  bottom-level  equation  can  be  expressed 
as 


where  J  is  the  circuit  Jacobian  matrix.  ecm(l)  0))  is  the  unit  vector  that  selects 
the  cosine  (sine)  part  of  the  first  harmonic  of  the  mth  cell  input  node  to  which  the 
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probe  is  attached.  Icmprobe  ( rmprobe )  is  the  cosine  (sine)  part  of  the  corresponding 
probe  current. 


The  top  level  solves  the  probe  equations  (Eq.  (5.15))  where  the  probe 
voltages  and  frequency  are  updated  until  convergence  is  reached  when  the 
currents  though  all  the  voltage  probes  go  to  zero. 


Ic  (yc  ...  yc  ys 

1,  probe  V  1,  probe  9  9  m,  probe  9  m,  probe  9 

Is  (yc  ...  Yc  Vs 

1,  probe  ^  1,  probe  9  9  m,  probe  9  m,  probe  9 


co)  =  0 


co)  =  0 


lc  (yc  ...  Yc  Vs 

m,  probe  '1,  probe  9  9  m,  probe  9  m,  probe  9 

Is  (yc  •••  Yc  Vs 

m,  probe  '1,  probe  9  9  m,  probe  9  m,  probe  9 


co)  =  0 


co)  =  0 


(5.15) 


The  Jacobian  for  this  set  of  equations  is: 
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(5.16) 


This  2 mx2m  Jacobian  can  be  obtained  by  a  simple  computation  using  information 
from  the  bottom- level  solution. 


The  frequency-domain  nonlinear  equations  for  the  circuit-probe 
combination,  G,  can  be  expressed  as: 


G(XA,a>)  =  0 


(5.17) 
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where  xA  is  the  vector  of  circuit  unknowns  and  the  probe  currents.  The  sensitivity 
with  respect  to  co  and  each  probe  voltage  can  then  be  computed  as: 


a  g 

dxA  _ 
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where  -  is  the  augmented  Jacobian  in  Eq.  (5.14).  Since  the  augmented 

dXA 

Jacobian  matrix  is  already  factored  in  the  bottom  level,  only  forward  and 
backward  substitutions  are  needed  to  solve  Eq.  (5.18).  All  entries  of  the  2 mx2m 
Jacobian  matrix  in  Eq.  (5.16)  can  be  easily  extracted  from  the  solution  of  Eq. 
(5.18). 


5. 5. 4  Computational  Effort  Analysis 

Since  multiple  probes  are  applied  to  the  circuit,  there  is  a  need  to  analyze  the 
computational  effort.  The  total  computational  effort  for  the  two-level  method  is 
dominated  by  the  bottom  level  due  to  two  reasons.  First,  the  bottom-level  circuit 
unknowns  are  much  larger  in  number  than  the  top-level  probe  unknowns.  Second, 
in  the  general  case,  the  number  of  bottom-level  iterations  is  several  times  larger 
than  the  number  of  top-level  iterations. 

From  Eq.  (5.14)  we  can  see  that  although  multiple  voltage  probes  have  been 
used,  the  bottom  level  matrix  size  increases  only  by  2(m-l)  (where  m  is  the 
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number  of  probes)  which  is  negligible  compared  to  the  bottom-level  matrix  size. 
Therefore,  the  increase  in  computational  effort  for  each  bottom-level  iteration  is 
minimal.  The  multiple-probe  method  is  almost  no  more  expensive  than  the 
traditional  single-probe  method  for  each  bottom-level  iteration. 

The  single-delay  cell  method  and  the  multiple-probe  method  have  been 
tested  on  a  wide  variety  of  single-ended  and  differential  ring  oscillators  to  show 
their  performance.  The  simulation  results  are  summarized  in  Chapter  7. 


5.6  Summary 

In  this  chapter  two  novel  ways  of  simulating  ring  oscillators  with  the 
harmonic  balance  method  are  described.  Both  of  the  approaches  exploit  the  phase 
relationships  of  a  ring  oscillator.  Details  of  the  numerical  formulation  and 
implementation  are  given.  The  proposed  single-delay  cell  method  assumes  that  all 
the  oscillator  cells  are  identical  and  cannot  be  used  when  large  mismatches 
between  cells  need  to  be  accounted  for.  The  multiple-probe  method  is  a  more 
general  method  which  can  be  used  instead. 


6.  FREQUENCY  DOMAIN  SIMULATION  OF  RELAXATION 

OSCILLATORS 
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The  operation  of  relaxation  oscillators  is  based  on  charge  and  discharge  of 
capacitors  with  a  constant  current  between  two  threshold  levels.  The  circuit 
waveforms  have  very  sharp  transitions.  Therefore  harmonic  truncation  produces 
large  aliasing  errors  which  makes  the  standard  harmonic  balance  method  not 
suitable  for  the  simulation  of  relaxation  oscillators. 


6.1  Relaxation  Oscillator  Frequency  Spectrum 

The  time-domain  waveform  as  well  as  the  frequency  domain  spectrum  for  a 
simple  CMOS  relaxation  oscillator  (Fig.  6.1)  are  shown  in  Fig.  6.2.  From  the 
figure  we  can  see  that  one  of  the  circuit  waveforms  has  a  very  sharp  transition  and 
the  signal  power  decreases  very  slowly  with  an  increase  in  the  number  of 
harmonics.  In  theory,  for  waveforms  with  sharp  transitions  the  Fourier 

coefficients  are  proportional  to  ,  where  n  is  the  number  of  harmonics.  So 

harmonic  truncation  will  introduce  significant  aliasing  error. 


Vdd 


Figure  6.1:  Schematic  of  a  CMOS  relaxation  oscillator. 
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Figure  6.2:  Time-domain  waveform  and  frequency  spectrum  for  a  relaxation 
oscillator. 


6.2  Aliasing  Error  for  Relaxation  Oscillator  Simulation 

The  relaxation  oscillator  in  Fig.  6.1  is  simulated  with  the  harmonic  balance 
method  to  show  the  influence  of  the  aliasing  error.  The  oscillation  frequency 
solution  from  the  harmonic  balance  method  using  a  different  number  of 
harmonics  is  shown  in  Fig.  6.3.  From  the  figure  we  can  see  that  with  an  increased 
number  of  harmonics  in  the  simulation,  the  solution  gets  closer  to  the  actual 
solution.  However,  even  with  more  than  100  harmonics,  it  is  still  off  by  a  factor  of 
2.  The  frequency  spectrum  of  the  circuit  waveform  is  shown  in  Fig.  6.4,  along 
with  the  result  from  a  transient  analysis.  The  discrepancy  between  the  results  from 
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the  harmonic  balance  method  and  the  transient  analysis  shows  the  extent  of 
aliasing  errors.  Fig.  6.5  gives  a  clear  picture  how  this  aliasing  error  influences  the 
time-domain  circuit  waveform  solution. 
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Figure  6.3:  Oscillation  frequency  solution  from  the  harmonic  balance  method 
using  a  different  number  of  harmonics. 


No.  of  Harmonics 


Figure  6.4:  Frequency  spectrum  from  transient  analysis  and  harmonic  balance 
method  with  128  harmonics. 
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Figure  6.5:  Time-domain  waveform  from  the  harmonic  balance  method  with  128 
harmonics. 


One  observation  is  that  to  faithfully  represent  the  circuit  waveforms  and 
effectively  control  the  aliasing  error,  harmonic  truncation  needs  to  take  into 
consideration  harmonics  whose  magnitude  is  greater  than  about  one-thousandth 
that  of  the  fundamental.  For  relaxation  oscillators  this  means  that  more  than  a 
thousand  harmonics  need  to  be  chosen  for  moderate  accuracy.  So  the  standard 
harmonic  balance  method  is  not  suitable  for  simulation  of  relaxation  oscillators. 
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7.  EXAMPLES  AND  RESULTS 


7.1  Frequency-Domain  Simulation  of  High-Q  Oscillators 

The  new  homotopy  algorithms  developed  in  Chapter  4  are  used  to  solve 
different  high-Q  oscillators  where  the  conventional  two-level  method  fails.  The 
results  from  the  homotopy-harmonic  balance  combined  method  are  verified  with 
the  results  from  long  and  expensive  transient  simulations.  The  developed  method 
is  also  used  to  simulate  moderate  Q  oscillators  and  then  the  computational 
performance  of  the  two-level  method  and  the  homotopy  method  is  compared. 


7.1.1  High-Q  Oscillator  Examples 

The  first  example  is  the  Pierce  oscillator  as  shown  in  Fig.  4.2.  The  voltage 
probe  is  attached  at  the  high  impedance  collector  node  of  the  BJT.  The  frequency 
from  a  pole-zero  analysis  is  used  as  the  initial  frequency  and  a  lmV  initial  value 
is  used  for  the  probe  voltage.  For  this  example,  the  whole  process  takes  9 
homotopy  steps  for  li.  The  total  number  of  iterations  (Jacobian  evaluations)  for 
Algorithm  1  and  Algorithm  2  is  942  and  304,  respectively.  The  oscillation 
frequency  is  found  to  be  9. 998422MHz.  The  solution  tracks  for  the  probe  voltage 
and  the  corresponding  frequency  using  Algorithm  1  and  Algorithm  2  are  the  same 
as  shown  in  Fig.  7.1.  From  the  figure  we  can  see  that  both  the  frequency  and  the 
probe  voltage  solution  tracks  rapidly  progress  towards  the  solution.  From  1  =0  to 
about  1=0.2,  the  frequency  changes  quickly,  and  the  tracking  takes  much  fewer 
steps  compared  with  that  of  Fig.  4.4.  This  is  a  clear  benefit  obtained  by 
decoupling  the  frequency  and  the  probe  voltage. 
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Figure  7.1:  Solution  traces  for  the  Pierce  oscillator  using  the  new  algorithm,  (a) 
Probe  voltage  as  a  function  of  X.  (b)  Oscillation  frequency  as  a  function  of  X. 

The  detailed  information  on  the  number  of  steps  at  Step  1  (Fig.  4.5  and  Fig. 
4.7)  for  each  A,i  using  Algorithm  1  and  Algorithm  2,  respectively,  is  summarized 
in  Table  7.1.  Algorithm  2  provides  a  speed-up  over  Algorithm  1  since  Newton’s 
method  instead  of  the  homotopy  method  is  used  to  solve  for  the  frequency.  This 
illustrates  that  Algorithm  2  is  efficient  and  it  is  used  for  all  subsequent  examples 
in  this  section. 
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Table  7.1:  Newton  iteration  information  for  each  7,i  step  for  the  Pierce  oscillator 
using  Algorithm  1  and  Algorithm  2.  AAy  is  the  Ay  increment  from  the  previous 
Ay  point. 


Xi  step 

A  A  i 

Algorithm  1 

Algorithm  2 

#  7.2  steps 

#  Newton  steps 

1 

0.0006 

22 

7 

2 

0.012 

31 

8 

3 

0.0047 

21 

7 

4 

0.10 

18 

6 

5 

0.20 

16 

6 

6 

0.27 

14 

5 

7 

0.32 

12 

6 

8 

0.36 

13 

4 

9 

-0.32 

9 

3 

To  verify  the  accuracy  of  the  homotopy  method,  a  long  and  expensive 
transient  analysis  with  tight  tolerances  is  performed  for  the  Pierce  oscillator.  The 
CPU  time  is  about  1000  times  that  of  the  homotopy  based  harmonic  balance 
method.  The  steady-state  output  voltage  in  both  the  time  and  frequency  domains 
is  shown  in  Fig.  7.2  along  with  the  results  from  a  transient  analysis.  From  the  plot 
we  can  see  that  the  results  from  the  homotopy  based  harmonic  balance  analysis 
and  transient  analysis  are  in  good  agreement. 

To  show  the  robustness  of  this  method  several  high-Q  oscillators  were 
simulated.  In  the  following  examples,  the  two-level  method  failed  because  an 
initial  guess  for  the  probe  voltage  could  not  be  determined  whereas  the  new 
homotopy  based  algorithm  converged.  Table  7.2  shows  the  starting  frequency/,,, 
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the  ratio  of  the  final  frequency  /  to  fa,  the  number  of  homotopy  steps,  the  total 
number  of  iterations,  and  the  CPU  times  on  a  SUN  UltralO  workstation.  Fig.  7.3 
and  Fig.  7.4  show  the  schematic  and  tracking  curves  for  the  Pierce2  oscillator. 
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Figure  7.2:  Results  from  transient  and  homotopy  based  harmonic  balance 
simulations  of  the  Pierce  oscillator. 


81 


Table  7.2:  Results  for  various  high-Q  oscillators. 


Circuit 

Q 

fo  (MHz) 

f/fo 

# 

Steps 

#Iter 

CPU  time 
(sec) 

Colpitts 

l.OOxlO3 

1.5915 

0.99987 

11 

208 

17.5 

TNT 

1.02xl03 

11.795 

0.84738 

15 

492 

26.3 

Pierce  (BJT) 

1.07xl04 

4.0811 

0.97036 

10 

241 

13.2 

Pierce  1  (MOS) 

3.69xl04 

7.9992 

0.99907 

5 

98 

8.2 

Clapp 

1.14xl05 

20.124 

0.99998 

12 

243 

19.2 

Pierce2  (MOS) 

2.82xl05 

1.1256 

0.99995 

26 

403 

141 

Figure  7.3:  Circuit  schematic  of  the  Pierce  oscillator  with  MOS  inverter. 
R1=1MQ,  R2=2.7Kn,  Cl=55pF,  C2=60pF,  Cp=25pF,  Cc=10fF,  Rc=60Q, 
Lc=2H. 
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Figure  7.4:  Solution  traces  for  the  high  Q  Pierce  oscillator  (Pierce2)  using  the  new 
algorithm,  (a)  Probe  voltage  as  a  function  of  X i .  (b)  Oscillation  frequency  as  a 
function  of  Xi. 


Oscillator  simulation  with  numerical  models  is  even  more  difficult  than  the 
analytical  models  as  discussed  in  [23].  The  new  algorithm  is  also  robust  in  this 
case.  Algorithm  2  has  been  implemented  in  the  coupled  device  and  circuit 
simulator,  CODECS  [45].  For  the  Pierce  1  oscillator  in  Table  7.2.  The  tracking 
takes  7  \  step  with  298  total  iterations.  The  solution  track  is  shown  in  Fig.  7.5. 
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Figure  7.5:  Solution  traces  for  the  Pierce  oscillator  (Piercel)  with  numerical 
models  using  the  new  algorithm,  (a)  Probe  voltage  as  a  function  of  (b) 
Oscillation  frequency  as  a  function  of  Xi. 


7.1.2  Moderate-Q  Oscillator  Examples 

The  proposed  homotopy  method  is  not  limited  to  high-Q  oscillators  only. 
Oscillators  with  moderate  Q  have  also  been  tested  [41].  The  results  are 
summarized  in  Table  7.3.  The  total  number  of  iterations  required  for  the  two- 
level  method  (for  the  initial  guess  and  for  solving  Eqs.  (3.8)  and  (3.9))  is  smaller 
than  that  for  the  homotopy  method  when  the  two-level  method  converges.  In 
these  situations,  the  homotopy  method  is  less  than  10  times  expensive  in  the  worst 
case,  and  for  most  circuits  it  is  only  2-4  times  slower.  However,  the  two-level 
method  does  not  converge  for  oscillator  circuits  with  numerical  models  whereas 
the  homotopy  method  does  (Table  7.3).  These  examples  illustrate  that  the  new 
homotopy-based  algorithm  is  robust  for  a  wide  range  of  circuits  and  models. 
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Table  7.3:  Comparison  of  the  two-level  method  and  homotopy  method  for 
moderate-Q  oscillator  circuits.  #  indicates  circuits  with  numerical  models,  and  x 
denotes  circuits  that  fail  to  converge. 


Circuit 

#  Iterations 

CPU  time  (sec) 

Two-level 

Homotopy 

Two- 

level 

Homotopy 

Eqs.  (5), 
(6) 

Colpitts  (BJT) 

44 

19 

197 

1.6 

5.2 

Colpitts  (MOS) 

52 

21 

208 

2.1 

6.3 

TNT 

113 

57 

310 

12.6 

23.4 

Wien 

64 

29 

274 

6.1 

14.5 

Sony 

56 

32 

170 

5.6 

12.1 

Phase  shift 

40 

17 

307 

23.9 

127.2 

Source  coupled 

58 

16 

358 

12.8 

60.1 

Emitter  coupled 
Colpitts 

43 

16 

256 

7.5 

36.7 

Cross  coupled 

45 

24 

398 

3.1 

20.8 

Colpitts(MOS)  # 

X 

X 

172 

X 

854.4 

Sony  # 

X 

X 

593 

X 

581.2 

7.2  Frequency-Domain  Simulation  of  Highly  Nonlinear  Ring  Oscillators 

In  this  section,  the  two  new  proposed  methods,  the  single-delay  cell  method 
and  the  multiple-probe  method  are  evaluated  on  a  variety  of  single-ended  and 
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differential  ring  oscillator  circuits.  The  convergence  performance  as  well  as  the 
accuracy  of  the  two  methods  is  demonstrated.  The  single-delay  cell  method 
provides  a  significant  speed-up  over  the  conventional  approach  as  shown.  The 
multiple-probe  method  is  efficient  and  converges  robustly  for  a  wide  variety  of 
ring  oscillator  circuits  with  non-identical  delay  cells. 


7.2.1  Convergence  Performance  and  Accuracy 

To  show  the  robustness  of  the  proposed  single-delay  cell  and  multiple-probe 
methods,  circuits  from  Table  5.2  have  been  simulated.  The  results  are  summarized 
in  Table  7.4.  From  the  results  we  can  see  that  the  two  proposed  methods  converge 
robustly  for  those  circuits  where  the  traditional  approach  fails.  The  two 
approaches  show  nearly  the  same  convergence  properties  and  the  convergence  is 
not  sensitive  to  the  initial  probe  voltage  and  the  oscillation  frequency  guess. 

The  frequency  initial  guess  is  estimated  from  the  time-delay  of  the  delay 
cells,  which  can  be  obtained  from  the  step-response  of  the  delay  cell.  The  probe 
voltage  initial  guess  is  simply  chosen  based  on  the  power  supply  voltage.  For  full¬ 
swing  delay  cells,  such  as  those  in  Circuit  1,  this  is  a  good  initial  guess  since  the 
output  swing  is  constrained  by  the  power  supply.  However,  for  partial-swing 
delay  cells  (Circuits  2,  3,  4)  the  circuit  structure  also  plays  a  role,  and  the  probe 
voltage  guess  depends  on  the  circuit.  The  new  single-delay  cell  and  multiple 
probe  methods  converge  robustly  for  ring  oscillator  circuits  with  partial-swing 
delay  cells  even  with  the  power  supply  voltage  used  as  the  initial  probe  voltage. 
This  characteristic  makes  the  application  of  the  harmonic  balance  method  to  ring 
oscillator  simulation  even  easier  since  the  only  initial  guess  that  needs  to  be 
specified  is  the  oscillation  frequency.  Furthermore,  as  shown  in  the  Table  7.4  this 
can  be  off  by  more  than  50%. 
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From  Table  7.4  we  can  see  that  the  convergence  of  the  ring  oscillator  top- 
level  probe  equation  is  not  sensitive  to  the  initial  guess  whenever  the  bottom- level 
iteration  converges.  This  is  consistent  with  the  low  Q  nature  of  ring  oscillator 
circuits.  However,  if  the  conventional  method  [7]  is  used  and  the  initial  guess  is 
too  far  from  the  real  solution,  the  bottom-level  equations  fail  to  converge  before 
the  frequency  and  probe  voltage  can  even  be  updated  at  the  top  level. 


Table  7.4:  Convergence  sensitivity  to  initial  guess  for  different  ring  oscillator 
circuits  using  the  new  single-delay  cell  and  multiple-probe  methods. 


Delay  cell 

Circuit  1 

Circuit  2 

Circuit  3 

Circuit  4 

|A/|//_ 

>50% 

>50% 

>50% 

>50% 

/vosc 

19.5% 

33.5% 

46.7% 

24.5% 

To  show  the  convergence  process  of  the  new  methods,  the  same  circuit 
waveform  (as  that  in  Fig.  5.2)  during  the  Newton  iterative  process  using  the 
single-delay  cell  and  multiple-probe  methods  is  shown  in  Fig.  7.6.  It  can  be  seen 
that  both  methods  constrain  the  circuit  waveform  while  requiring  no  damping  or 
any  thresholds.  Also  the  convergence  is  achieved  in  only  a  few  iterations. 
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(a) 


(b) 

Figure  7.6:  Circuit  waveform  during  the  Newton  iterative  process,  (a)  Single¬ 
delay  cell  method,  (b)  Multiple-probe  method. 


The  steady-state  output  voltage  in  both  the  time  and  frequency  domains  for 
Circuit  1  is  shown  in  Fig.  7.7  along  with  the  results  from  a  transient  analysis  to 
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show  the  accuracy  of  the  new  methods.  From  the  plots  we  can  see  that  the  results 
from  the  single-delay  cell  and  multiple-probe  harmonic  balance  methods  are  in 
good  agreement  with  those  from  the  transient  analysis. 
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Figure  7.7:  Results  from  transient,  single-delay  cell  harmonic  balance,  and 
multiple-probe  harmonic  balance  simulations  of  a  9-stage  ring  oscillator. 


7.2.2  Computational  Performance  of  the  Single-delay  Cell  Method 


The  new  single-delay  cell  method  is  very  efficient  for  ring  oscillator 
simulation.  Several  single-ended  and  differential  ring  oscillators  have  been 
simulated  using  the  conventional  approach  and  the  new  single-delay  cell 
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approach.  The  inverters  are  implemented  by  the  current- starved  delay  cell  (Fig. 
5.1(b))  and  the  cross-coupled  load  delay  cell  (Fig.  5.1(d)),  respectively.  From  the 
simulation  results  in  Table  7.5  we  can  see  that  the  single-delay  cell  method 
provides  a  significant  speed-up  in  simulations  even  for  a  3  stage  differential  ring 
oscillator.  The  number  of  iterations  for  the  single-delay-cell  method  is  also  less 
than  the  complete  oscillator  simulation  due  to  a  smaller  number  of  nonlinear 
devices. 


Table  7.5:  Comparison  of  the  conventional  approach  with  the  new  single-delay 
cell  approach  for  ring  oscillators. 


Circuit 

Method 

#Unknowns 

Matrix  size 

# 

Iterations 

Memory 

(M) 

CPU 

time(min) 

1 

Conventional 

19 

933x933 

201 

20.7 

20.25 

Delay  cell 

7 

345x345 

97 

8.84 

0.70 

2 

Conventional 

69 

1727x1727 

355 

18.2 

39.63 

Delay  cell 

31 

777x777 

131 

11.3 

2.55 

Circuit  1:  15-stage  single-ended  ring  oscillator  with  current  starved  delay  cell. 
Circuit  2:  3-stage  differential  ring  oscillator  with  cross-coupled  load  delay  cell. 

One  key  benefit  of  the  single-delay  cell  method  is  that  ring  oscillators  with  a 
different  number  of  stages  can  be  simply  simulated  by  using  the  single-delay  cell 
equivalent  circuit  with  the  value  of  the  phase  delay  0  adjusted  according  to  Eq. 
(5.2).  The  simulation  results  are  presented  in  Tables  7.6  and  7.7  for  a  single-ended 
ring  oscillator  and  a  differential  ring  oscillator,  respectively.  The  current- starved 
delay  cell  (Fig.  5.1(b))  and  the  Maneatis  delay  cell  [42]  (Fig.  7.9(a))  are  used. 
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From  the  tables,  we  see  that  the  more  the  number  of  stages,  the  higher  the 
nonlinearity,  and  therefore,  the  more  the  CPU  time  required  for  convergence. 
However,  for  the  single-delay  cell  method,  the  CPU  time  only  increases  by  a 
factor  of  2-3  whereas  for  the  conventional  method,  the  simulation  time  increases 
by  a  factor  of  40  for  the  21 -stage  single-ended  ring  oscillator  and  the  convergence 
is  very  sensitive  to  the  initial  guess.  For  differential  ring  oscillators  with  more 
than  3  stages  the  conventional  method  failed  to  reach  convergence.  The  slight 
increase  in  the  CPU  time  for  the  single-delay  cell  method  with  an  increased 
number  of  stages,  is  due  to  the  circuit  operation  becoming  more  and  more 
nonlinear. 


Table  7.6:  Comparison  of  simulation  results  for  single-ended  ring  oscillators  with 
different  number  of  stages. 


#Stages 

5 

7 

11 

15 

17 

19 

21 

CPU  time 
(min) 

Single-delay 
Cell  Method 

0.15 

0.18 

0.20 

0.23 

0.25 

0.28 

0.42 

Conventional 

Method 

0.62 

2.60 

4.12 

6.55 

10.92 

17.88 

25.35 

Memory 

(M) 

Single-delay 
Cell  Method 

7.2 

7.2 

7.2 

7.2 

7.2 

7.2 

7.2 

Conventional 

Method 

8.9 

9.7 

11.2 

12.6 

13.2 

14.1 

15.3 
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Table  7.7:  Comparison  of  simulation  results  for  differential  ring  oscillators  with 
different  number  of  stages,  x  denotes  no  convergence  under  exact  initial  guess  for 
the  oscillation  frequency  and  probe  voltage. 


#Stages 

3 

5 

7 

9 

CPU  time  (min) 

Single-delay  Cell 
Method 

2.31 

3.22 

4.00 

4.85 

Conventional 

Method 

34.63 

X 

X 

X 

Memory  (M) 

Single-delay  Cell 
Method 

10.9 

10.9 

10.9 

10.9 

Conventional 

Method 

17.2 

23.5 

35.7 

53.4 

The  singe  delay  cell  method  is  also  effective  for  differential  ring  oscillators 
with  an  even  number  of  stages  since  the  phase  relationship  for  each  delay  cell  still 
holds.  A  4-stage  differential  ring  oscillator  (Fig.  7.8)  is  simulated  with  three  kinds 
of  delay  cells:  Maneatis  delay  cell  [42],  Lee/Kim  delay  cell  [43]  and  current- 
controlled  delay  cell  [44]  (Fig.  7.9).  All  circuits  converge  robustly  and  the 
performance  is  summarized  in  Table.  7.8.  There  is  also  a  good  agreement  between 
the  results  from  the  harmonic  balance  method  and  transient  analysis  as  shown  in 
Figure  7.8  for  the  current-controlled  delay  cell. 


Figure  7.8:  4-stage  differential  ring  oscillator. 
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Table  7.8:  Simulation  performance  comparison  for  the  four-stage  differential  ring 
oscillators  with  different  delay  cells  using  the  single-delay  cell  approach. 


Delay  cells 
(Fig.  7.9) 

#Unknowns 

Matrix  size 

#  Iterations 

CPU  time(min) 

(a) 

29 

727x727 

158 

3.17 

(b) 

22 

552x552 

201 

2.37 

(c) 

32 

802x802 

85 

1.67 

0  5  10  15  20  25  30  35  40 

No. of  Harmonics 


Figure  7.10:  Simulation  results  for  the  4-stage  differential  ring  oscillator  (current- 
controlled  delay  cell). 
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7.2.3  Computational  Performance  of  the  Multiple-probe  Method 

For  ring  oscillator  circuits  that  converge  with  the  conventional  single  probe 
method  using  a  close  initial  guess,  the  multiple-probe  method  is  more  efficient 
even  when  more  probes  are  applied  in  the  circuit.  The  simulation  results  are 
shown  in  Table. 7. 9.  In  this  case  the  matrix  size  and  memory  requirements  for  the 
two  methods  are  similar.  However,  the  multiple-probe  method  uses  a  more 
effective  way  to  force  the  circuit,  thereby  requiring  fewer  iterations  for 
convergence. 


Table  7.9:  Comparison  of  the  multiple-probe  method  with  the  single-probe 
method  for  ring  oscillator  simulation. 


Cir 

# 

Stages 

Delay  cell 

#  Iterations 

(#  Top  level  iterations) 

CPU  time  Multiple- 
probe/Conventional 
Method 

Multiple-Probe 

Method 

Conventional 

Method 

1 

7 

Fig.  5.1(a) 

58  (6) 

104  (5) 

0.53 

2 

9 

Fig.  5.1(b) 

48  (4) 

75  (4) 

0.63 

3 

3 

Fig.  5.1(d) 

39  (5) 

56(5) 

0.67 

4 

3 

Fig.  7.9(a) 

64  (4) 

90  (4) 

0.71 

5 

4 

Fig.  7.9(b) 

59(6) 

128  (7) 

0.46 
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7.2.4  Multiple-probe  Method  for  Non-ideal  Delay  Cells 

The  single-delay  cell  method  is  very  efficient  as  shown  in  Section  7.2.2. 
However,  it  cannot  handle  ring  oscillators  with  non-identical  delay  cells.  The 
multiple-probe  method  simulates  the  full  ring  oscillator  circuit.  No  assumptions 
about  identical  delay  cells  are  made,  so  it  can  be  used  to  simulate  non-identical 
delay  cells.  Next  the  effectiveness  of  the  multiple-probe  method  for  simulating 
ring  oscillators  with  mismatched  or  non-identical  delay  cells  is  demonstrated.  The 
results  for  several  circuits  are  summarized  in  Table  7.10.  The  first  circuit  is 
connected  to  an  output  buffer  and  the  second  circuit  was  implemented  using 
mismatches  within  the  delay  cells.  For  simplicity  we  assume  the  mismatches  only 
exist  between  the  input  NMOSFET  of  each  stage.  From  the  table  we  can  see  that 
the  proposed  method  converges  robustly  for  these  cases.  The  number  of  iterations 
with  identical  delay  cell  cases  and  mismatch/loading  cases  are  quite  similar. 


Table  7.10:  Comparison  of  ring  oscillator  simulation  with  identical  delay  cells  and 
mismatched  cells  using  the  multiple-probe  method. 


Delay  cell 

Fig.  5.1(a) 

Fig.  7.9(c) 

#  Stages 

5 

4 

Oscillator  circuit 

Without 
output  buffer 

With  output 
buffer 

Identical 

cells 

Mismatched 

cells 

Oscillation  frequency 
(MHz) 

237.56 

216.74 

116.92 

121.82 

#Iter 

Top  level 

8 

7 

5 

5 

Bottom  level 

53 

47 

70 

76 
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Though  the  initial  guesses  for  multiple  probes  are  based  on  the  assumption 
of  identical  delay  cells,  the  proposed  method  is  able  to  handle  ring  oscillators  with 
large  mismatches  in  the  delay  cells.  To  demonstrate  the  robustness  of  the 
multiple-probe  method,  an  extreme  case  is  used.  A  nine- stage  stepped  buffer  is 
connected  as  a  ring  oscillator  configuration  as  shown  in  Fig.  7.11.  The  multiple- 
probe  method  converges  and  requires  only  106  bottom  level  and  11  top  level 
iterations  while  the  single-probe  method  fails. 


Figure  7.11:  9-stage  stepped-buffer  connected  as  a  ring  oscillator.  Each  stage  is 
sized  a  factor  of  e  larger  than  the  previous  stage. 


7.2.5  Simulation  of  Ring  Oscillators  with  Numerical  Models 

In  this  section,  more  accurate  numerical  models  are  used  for  ring  oscillator 
simulation.  The  new  algorithms  for  robust  ring  oscillator  simulation  are 
implemented  in  the  coupled  device  and  circuit  simulator  CODECS  [24,  45].  A 
two-dimensional  (2D)  numerical  MOSFET  model  with  31x19  mesh  points  is  used. 
While  the  single-probe  method  fails  for  these  circuits  with  numerical  models,  the 
single-delay  cell  and  multiple-probe  methods  deliver  reliable  convergence.  The 
simulation  results  are  the  same  for  the  two  methods.  Since  circuit  simulation  with 
numerical  models  is  very  time-consuming,  the  single-delay  cell  method  is 
preferred  in  these  cases. 
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Table  7.11:  Ring  oscillator  simulation  results  with  2D  numerical  MOSFET 
models. 


Delay  cells 

Single-delay  Cell  Method 

Multiple-probe  Method 

#  Iterations 
(#  Top  level 
iterations) 

CPU 

Time 

(min) 

Memory 

(M) 

#  Iterations 
(#  Top  level 
iterations) 

CPU 

Time 

(min) 

Memory 

(M) 

Fig.  5.1  (d) 

270(5) 

54.4 

14.3 

237  (5) 

164.9 

35.5 

Fig.  7.9  (a) 

98(5) 

18.1 

15.4 

101  (5) 

54.60 

31.9 

Fig.  7.9  (b) 

136(5) 

17.4 

13.9 

93  (4) 

52.37 

35.4 

Fig.  7.9  (c) 

121(4) 

91.2 

20.8 

129  (6) 

224.4 

44.8 

The  influence  of  substrate  noise  on  the  output  spectrum  of  the  ring  oscillator 
has  also  been  studied  with  more  accurate  numerical  models.  The  set-up  for  the 
simulation  is  as  shown  in  Fig.  7.12.  There  are  two  tones  in  the  circuit.  One  is  the 
unknown  oscillator  frequency,  the  other  is  the  substrate  noise  frequency.  The 
problem  thus  becomes  a  mixer  plus  oscillator  problem.  In  this  particular  example 
52  harmonics  were  chosen  by  the  box  truncation  method  [2].  The  simulation 
results  with  and  without  substrate  noise  are  shown  in  Fig.  7.13.  From  the  figure 
we  can  clearly  see  the  effect  of  the  substrate  noise  on  the  oscillator  output 
spectrum. 


Magnitude  (dB) 


Figure  7.12:  Ring  oscillator  simulation  with  substrate  noise.  fosc=0.414GHz, 
fsub=0.8GHz,  Vsub=10mV. 


Figure  7.13:  Simulation  results  for  a  3-stage  ring  oscillator  with  substrate  noise. 
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7.3  Summary 

In  this  chapter,  oscillator  examples  are  used  to  show  the  performance  of  the 
proposed  homotopy-based  harmonic  balance  method,  single-delay  cell  method 
and  the  multiple  probe  method. 

The  homotopy-based  harmonic  balance  method  is  robust  and  efficient  for 
the  simulation  of  high-Q  oscillators  where  the  conventional  two-level  method 
failed.  This  method  can  also  be  used  for  the  simulation  of  moderate  Q  oscillators 
with  acceptable  computational  cost.  While  the  two-level  method  fails  for  some 
oscillator  circuits  with  numerical  models,  the  homotopy-based  harmonic  balance 
method  delivers  reliable  convergence. 

The  single-delay  cell  and  multiple-probe  methods  solve  the  convergence 
problems  in  the  simulation  of  ring  oscillators.  Both  methods  are  robust  and  not 
sensitive  to  the  initial  guesses.  The  single-delay  cell  method  provides  a  significant 
speed-up  over  the  multiple-probe  method. 
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8.  CONCLUSIONS 


8.1  Summary  of  Contributions 

This  work  presents  the  first  detailed  study  of  frequency-domain  simulation 
of  oscillators  using  the  harmonic  balance  method.  The  convergence  problems 
encountered  in  different  kinds  of  “hard-to-solve”  oscillator  circuits  have  been 
studied.  New  approaches  for  robust  and  efficient  harmonic  balance  simulation  of 
high-Q  oscillators  and  highly  nonlinear  ring  oscillators  have  been  developed.  The 
application  of  the  frequency-domain  harmonic  balance  method  to  relaxation 
oscillator  simulation  was  also  explored. 

A  globally  convergent  homotopy  method  is  applied  for  robust  high-Q 
oscillator  simulation  in  the  frequency  domain.  The  method  uses  the  idea  of  the 
voltage  probe  while  requiring  no  initial  guess  for  the  probe  voltage.  Different 
embedding  techniques  have  been  evaluated  and  an  efficient  and  robust  algorithm 
has  been  demonstrated.  A  variety  of  high-Q  oscillators  as  well  as  moderate-Q 
oscillators  were  simulated  to  show  the  validity  of  the  new  method.  The  developed 
algorithm  can  be  easily  applied  to  existing  harmonic  balance  software  by  linking 
with  the  HOMPACK  package. 

Two  novel  ways  of  simulating  ring  oscillators  with  the  harmonic  balance 
method  were  also  proposed.  The  phase  relationships  of  a  ring  oscillator  is 
exploited  in  both  approaches.  In  the  single-delay  cell  method,  a  single-delay  cell 
equivalent  circuit  is  developed  for  ring  oscillator  simulation.  It  is  shown  that  the 
phase  relationship  can  be  easily  and  efficiently  handled  in  the  frequency  domain. 
Details  of  the  numerical  formulation  and  implementation  are  given.  With  a 
reduced  number  of  nodes  and  nonlinear  devices  involved  in  the  simulation,  the 
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single-delay  cell  method  results  in  a  significant  improvement  in  the  simulation 
time  and  convergence.  This  approach  can  also  be  applied  to  a  flat  netlist  with  the 
aid  of  subcircuit  identification  algorithms  [46].  The  proposed  method  assumes 
that  all  the  oscillator  cells  are  identical  and  cannot  be  used  when  large  mismatches 
between  cells  need  to  be  accounted  for. 

The  proposed  multiple-probe  method  can  handle  ring  oscillator  circuits  with 
large  mismatches  in  the  delay  cells.  In  the  multiple-probe  method,  a  voltage  probe 
is  applied  to  the  input  of  the  each  delay  cell.  Although  more  than  one  probe  is 
used,  the  new  method  requires  no  additional  initial  guesses  compared  with  the 
conventional  single-probe  method.  The  use  of  multiple  probes  provides  a  method 
that  is  robust  for  ring  oscillator  simulation.  Furthermore,  this  method  is  also 
efficient  compared  with  the  conventional  method. 

A  variety  of  ring  oscillator  circuits  have  been  simulated  showing  the  validity 
of  the  two  proposed  approaches.  If  the  delay  cells  in  the  ring  oscillators  are 
identical,  the  single-delay  cell  method  is  recommended  for  use,  since  it  provides  a 
significant  speed  up  over  the  conventional  method  and  the  multiple-probe  method. 

The  standard  harmonic  balance  method  is  not  suitable  for  simulation  of  the 
relaxation  type  oscillators  due  to  accuracy  problems  arising  from  aliasing  errors. 


8.2  Future  Work 

For  some  ring  oscillators  with  complex  delay  cells,  the  frequency-domain 
harmonic  balance  simulation  results  in  a  large  system  matrix.  Direct  factorization 
of  such  a  large  matrix  is  time  consuming  and  Krylov  subspace  methods  should  be 
used.  However,  since  the  ring  oscillator  is  quite  nonlinear,  a  conventional  block 
diagonal  preconditioner  is  not  a  good  choice.  More  work  needs  to  be  done  to  find 
a  suitable  preconditioner. 
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In  the  single-delay  cell  method,  the  equivalent  circuit  used  for  simulation  is 
based  on  identical  delay  cells.  Extensions  of  this  method  to  non-ideal  delay  cells 
needs  further  investigation.  A  sensitivity-like  analysis  can  be  used  in  conjunction 
with  the  single-delay  cell  method  for  small  mismatches.  The  algorithmic  approach 
needs  to  be  developed  and  evaluated  for  both  accuracy  and  performance. 

It  has  been  shown  that  the  standard  harmonic  balance  method  is  not 
suitable  for  simulation  of  relaxation  oscillators.  Some  other  approaches  have  been 
proposed  to  handle  circuits  with  sharp  transitions,  such  as  the  time-mapped 
harmonic  balance  method  [47-48],  the  wavelet  method  [49-51]  and  multi-interval 
Chebyshev  method  [52].  However,  it  is  not  yet  clear  whether  extensions  of  these 
methods  to  the  simulation  of  relaxation  oscillators  is  possible.  The  combination  of 
frequency-domain  and  time-domain  methods  may  also  be  a  solution  and  needs 
further  investigation. 
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